Library Maths Release v2.1 -A library of classes and functions for mathematical objects.

1aTabOrder

Library Maths Tab Order Version: 2.1
------------------------------
This file should not be included in the Codea project.
#ChangeLog
#Main
#Matrix
#Complex
#Quaternion
#Vec3
#Vector
#ccConfig

ccConfig.lua

--[[
###########################################
##Codea Community Project Config Settings##
###########################################

##You can use # to comment out a line
##Use 1 for true and 0 for false


###########################################
#       Add project info below            #
#==========================================
ProjectName: Place Project Name Here
Version: Alpha 1.0
Comments:
Author:
##License Info: http://choosealicense.com
##Supported Licneses: MIT, GPL, Apache, NoLicense
License: MIT
#==========================================


###########################################
#                Settings                 #
[Settings]=================================
##Codea Community Configuration settings
##Format: Setting state

Button 1
NotifyCCUpdate 1
ResetUserOption 0
AddHeaderInfo 1

[/Settings]================================



###########################################
#              Screenshots                #
[Screenshots]==============================
##Screenshots from your project.
##Format: url
##Example: http://www.dropbox.com/screenshot.jpg


[/Screenshots]=============================



###########################################
#                   Video                 #
[Video]====================================
##Link to a YouTube.com video.
##Format: url
##Example: http://www.youtube.com/videolink


[/Video]===================================



###########################################
#              Dependencies               #
[Dependencies]=============================
##Include the names of any dependencies here
##Format: Dependency
##Example: Codea Community


[/Dependencies]============================



############################################
#                   Tabs                   #
[Tabs]======================================
##Select which tabs are to be uploaded.
##Keyword 'not' excludes a tab or tabs. Keyword 'add' includes a tab or tabs.
##not * will exclude all tabs to allow for individual selection
##not *tab1 will look for any tab with tab1 on the end.
##not tab1* will look for any tab with tab1 at the beginning.
##Format: (add/not)(*)tab(*)
##Example: not Main --this will remove main.
##Example: not *tab1 --this will remove any tab with "tab1" on the end.
##Example: add Main --This will add Main back in.



[/Tabs]=====================================



#############################################
#                  Assets                   #
[Assets]=====================================
##Directory, path and url info for any assets besides the standard Codea assets.
##Format: Folder:sprite URL
##Example: Documents:sprite1 http://www.somewebsite.com/img.jpg


[/Assets]====================================
--]]

ChangeLog.lua

--[[
ChangeLog
=========

v2.1 Fixed bug in Quaternion "tomatrix" function.

v2.0 Split into separate libraries: "Library Maths" contains the
     mathematical modules.

v1.0 Converted to use toadkick's cmodule for importing 
     and Briarfox's AutoGist for versionning.
     It needs toadkick's code from 
        https://gist.github.com/apendley/5411561
     tested with version 0.0.8
     To use without AutoGist, comment out the lines in setup.

--]]

Complex.lua

--[==[
-- Complex numbers class
-- 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 complex numbers.
--]]

local Complex = class()

Complex.symbol = readLocalData("Complex Symbol","i")
local radeg = readLocalData("Complex Angle","rad")
if radeg == "rad" then
    Complex.angle = 1
    Complex.angsym = "π"
else
    Complex.angle = 180
    Complex.angsym = "°"
end
Complex.precision = readLocalData("Complex Precision",2)


--[[
A complex number can either be specified by giving the two
coordinates as real numbers or by giving a vector.
--]]


function Complex:init(...)
    -- you can accept and set parameters here
    if arg.n == 2 then
        -- two numbers
        self.z = vec2(arg[1],arg[2])
    elseif arg.n == 1 then
        -- vector
        self.z = arg[1]
    else
        print("Incorrect number of arguments to Complex")
    end
end

function Complex:clone(z)
    return Complex(z:real(),z:imaginary())
end

--[[
Test if we are zero.
--]]

function Complex:is_zero()
    -- are we the zero vector
    if self.z ~= vec2(0,0) then
        return false
    end
    return true
end

--[[
Test if we are real.
--]]

function Complex:is_real()
    -- are we the zero vector
    if self.z.y ~= 0 then
        return false
    end
    return true
end

--[[
Test if the real part is zero.
--]]

function Complex:is_imaginary()
    -- are we the zero vector
    if self.z.x ~= 0 then
        return false
    end
    return true
end

--[[
Test for equality.
--]]

function Complex:is_eq(q)
    if self.z ~= q.z then
        return false
    end
    return true
end

--[[
Defines the "==" shortcut.
--]]

function Complex:__eq(q)
    return self:is_eq(q)
end

--[[
The inner product of two complex numbers.
Why did I program this?
--]]

function Complex:dot(q)
    return self.z:dot(q.z)
end

--[[
Length of a complex number.
--]]

function Complex:len()
    return self.z:len()
end

--[[
Often enough to know the length squared, which is quicker.
--]]

function Complex:lensq()
    return self.z:lenSqr()
end

--[[
Distance between two complex numbers.
--]]

function Complex:dist(w)
    return self.z:dist(w.z)
end

--[[
Often enough to know the distance squared, which is quicker.
--]]

function Complex:distSqr(w)
    return self.z:distSqr(w.z)
end

--[[
Normalise a complex number to have length 1, if possible.
--]]

function Complex:normalise()
    local l
    if self:is_zero() then
        print("Unable to normalise a zero-length complex number")
        return false
    end
    l = 1/self:len()
    return self:scale(l)
end

--[[
Scale the complex number.
--]]

function Complex:scale(l)
    return Complex(l * self.z)
end

--[[
Add two complex numbers.  Or add a real number to a complex one.
--]]

function Complex:add(q)
    if type(q) == "number" then
        return Complex(self.z + vec2(q,0))
    else
        return Complex(self.z + q.z)
    end
end

--[[
q + p
--]]

function Complex:__add(q)
    if type(self) == "number" then
        return q:add(self)
    else
        return self:add(q)
    end
end

--[[
Subtraction
--]]

function Complex:subtract(q)
    if type(q) == "number" then
        return Complex(self.z - vec2(q,0))
    else
        return Complex(self.z - q.z)
    end
end

--[[
q - p
--]]

function Complex:__sub(q)
    if type(self) == "number" then
        return q:subtract(self):scale(-1)
    else
        return self:subtract(q)
    end
end

--[[
Negation (-q)
--]]

function Complex:__unm()
    return self:scale(-1)
end

--[[
Length (#q)
--]]

function Complex:__len()
    return self:len()
end

--[[
Multiplication.
--]]

function Complex:multiply(q)
    local a,b,c,d
    a = self.z.x
    b = self.z.y
    c = q.z.x
    d = q.z.y
    return Complex(a*c - b*d,a*d + b*c)
end

--[[
q * p
--]]

function Complex:__mul(q)
    if type(q) == "number" then
        return self:scale(q)
    elseif type(self) == "number" then
        return q:scale(self)
    elseif type(q) == "table" then
        if q:is_a(Complex) then
                return self:multiply(q)
        end
    end
end

--[[
Conjugation.
--]]

function Complex:conjugate()
    return Complex(self.z.x, - self.z.y)
end

function Complex:co()
    return self:conjugate()
end

--[[
Reciprocal: 1/q
--]]

function Complex:reciprocal()
    if self:is_zero() then
        print("Cannot reciprocate a zero complex number")
        return false
    end
    local q = self:conjugate()
    local l = self:lensq()
    q = q:scale(1/l)
    return q
end

--[[
Real powers.
--]]

function Complex:repower(n,k)
    local r = self.z:len()
    local t = -self.z:angleBetween(vec2(1,0))
    k = k or 0
    r = math.pow(r,n)
    t = (t + k * 2 * math.pi) *n
    return Complex:polar(r,t)
end

--[[
Complex powers.
--]]

function Complex:power(w,k)
    if type(w) == "number" then
        return self:repower(w,k)
    end
    if self:is_zero() then
        print("Taking powers of 0 is somewhat dubious")
        return false
    end
    local r = self.z:len()
    local t = -self.z:angleBetween(vec2(1,0))
    local u = w.z.x
    local v = w.z.y
    k = k or 0
    local nr = math.pow(r,u) * math.exp(-v*t)
    local nt = (t + k * 2 * math.pi) * u + math.log(r) * v
    return Complex:polar(nr,nt)
end

--[[
q^n

This is overloaded so that a non-(complex) number exponent returns
the conjugate.  This means that one can write things like q^"" to
get the conjugate of a complex number.
--]]

function Complex:__pow(n)
    if type(n) == "number" then
        return self:repower(n)
    elseif type(n) == "table" and n:is_a(Complex) then
        return self:power(n)
    else
        return self:conjugate()
    end
end

--[[
Division: q/p
--]]

function Complex:__div(q)
    if type(q) == "number" then
        return self:scale(1/q)
    elseif type(self) == "number" then
        return q:scale(1/self):reciprocal()
    elseif type(q) == "table" then
        if q:is_a(Complex) then
                return self:multiply(q:reciprocal())
        end
    end
end

--[[
Returns the real part.
--]]

function Complex:real()
    return self.z.x
end

--[[
Returns the imaginary part.
--]]

function Complex:imaginary()
    return self.z.y
end

--[[
Represents a complex number as a string.
--]]

Complex.precision = 2

function Complex:__concat(v)
    if type(v) == "table" 
        and v:is_a(Complex) then
            return self .. v:tostring()
        else
            return self:tostring() .. v
        end
end

function Complex:__tostring()
    return self:tostring()
end

function Complex:tostring()
    local s
    local x = math.floor(
        self.z.x * 10^Complex.precision +.5
        )/10^Complex.precision
    local y = math.floor(
        self.z.y * 10^Complex.precision +.5
        )/10^Complex.precision
    if x ~= 0 then
        s = x
    end
    if y ~= 0 then
        if s then 
                if y > 0 then
                    if y == 1 then
                        s = s .. " + " .. Complex.symbol
                    else
                        s = s .. " + " .. y .. Complex.symbol
                    end
                else
                    if y == -1 then
                        s = s .. " - " .. Complex.symbol
                    else
                        s = s .. " - " .. (-y) .. Complex.symbol
                    end
                end
        else
                if y == 1 then
                    s = Complex.symbol
                elseif y == - 1 then
                    s = "-" .. Complex.symbol
                else
                    s = y .. Complex.symbol
                end
        end
    end
    if s then
        return s
    else
        return "0"
    end
end

function Complex:topolarstring()
    local t = math.floor(Complex.angle *
        self:arg() * 10^Complex.precision/math.pi +.5
        )/10^Complex.precision
    local r = math.floor(
        self:len() * 10^Complex.precision +.5
        )/10^Complex.precision
    local s = ""
    --[[
    -- this is exponential notation
    if t == 0 then
        s = r
    else
        if r ~= 1 then
            s = r
        end
        s = s .. "e^(" .. t .."i)"
    end
    --]]
    s = "(" .. r .. "," .. t .. Complex.angsym .. ")"
    return s
end

function Complex:arg()
    return -self.z:angleBetween(vec2(1,0))
end

function Complex:polar(r,t)
    return Complex(r * math.cos(t), r * math.sin(t))
end

--[[
The unit complex number.
--]]

function Complex.unit()
    return Complex(1,0)
end
function Complex.zero()
    return Complex(0,0)
end
function Complex.i()
    return Complex(0,-1)
end

--[[
Overload the maths functions
--]]

local __math = {}
for _,v in ipairs({
    "abs",
    "pow",
    "sqrt",
    "exp",
    "log",
    "sin",
    "cos",
    "sinh",
    "cosh",
    "tan",
    "tanh"
    }) do
    __math[v] = math[v]
end

function math.abs(n)
    if type(n) == "number" then return __math.abs(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            return n:len()
    end
    print("Cannot take the length of " .. n)
end

function math.pow(n,e)
    if type(n) == "number" then
        if type(e) == "number" then
             return __math.pow(n,e)
        elseif type(e) == "table"
        and e.is_a
        and e:is_a(Complex)
        then
            local w = Complex(n,0)
            return w:repower(e,0)
        end
    end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            return n:power(e,0)
    end
    print("Cannot take the power of " .. n .. " by " .. e)
end

function math.sqrt(n)
    if type(n) == "number" then return __math.sqrt(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            return n:repower(.5,0)
    end
    print("Cannot take the square root of " .. n)
end

function math.exp(n)
    if type(n) == "number" then return __math.exp(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            return Complex:polar(math.exp(n:real()),n:imaginary())
    end
    print("Cannot exponentiate " .. n)
end

--[[
cos(x+iy) = cos(x) cos(iy) - sin(x) sin(iy)
          = cos(x) cosh(y) - i sin(x) sinh(y)
--]]

function math.cos(n)
    if type(n) == "number" then return __math.cos(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            local x = n:real()
            local y = n:imaginary()
            return Complex(__math.cos(x)*__math.cosh(y),-__math.sin(x)*__math.sinh(y))
    end
    print("Cannot take the cosine of " .. n)
end

--[[
sin(x+iy) = sin(x) cos(iy) + cos(x) sin(iy)
          = sin(x) cosh(y) + i cos(x) sinh(y)
--]]

function math.sin(n)
    if type(n) == "number" then return __math.sin(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            local x = n:real()
            local y = n:imaginary()
            return Complex(__math.sin(x)*__math.cosh(y), __math.cos(x)*__math.sinh(y))
    end
    print("Cannot take the sine of " .. n)
end

--[[
cosh(x+iy) = cosh(x) cosh(iy) + sinh(x) sinh(iy)
           = cosh(x) cos(y) + i sinh(x) sin(y)
--]]

function math.cosh(n)
    if type(n) == "number" then return __math.cosh(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            local x = n:real()
            local y = n:imaginary()
            return Complex(__math.cosh(x)*__math.cos(y), __math.sinh(x)*__math.sin(y))
    end
    print("Cannot take the hyperbolic cosine of " .. n)
end

--[[
sinh(x+iy) = sinh(x) cosh(iy) + cosh(x) sinh(iy)
          = sinh(x) cos(y) + i cosh(x) sin(y)
--]]

function math.sinh(n)
    if type(n) == "number" then return __math.sinh(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            local x = n:real()
            local y = n:imaginary()
            return Complex(__math.sinh(x)*__math.cos(y), __math.cosh(x)*__math.sin(y))
    end
    print("Cannot take the hyperbolic sine of " .. n)
end

--[[
tan(x+iy) = (sin(x) cos(x) + i sinh(y) cosh(y))
            /(cos^2(x) cosh^2(y) + sin^2(x) sinh^2(y))
--]]

function math.tan(n)
    if type(n) == "number" then return __math.tan(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            local x = n:real()
            local y = n:imaginary()
            local cx = __math.cos(x)
            local sx = __math.sin(x)
            local chy = __math.cosh(y)
            local shy = __math.sinh(y)
            local d = cx^2 * chy^2 + sx^2 * shy^2
            if d == 0 then
                return false
            end
            return Complex(sx*cx/d,shy*chy/d)
    end
    print("Cannot take the tangent of " .. n)
end

--[[
tanh(x+iy) = i tan(y - ix)
           = (sin(x) cos(x) + i sinh(y) cosh(y))
            /(cos^2(x) cosh^2(y) + sin^2(x) sinh^2(y))
           = (sinh(x) cosh(x) + i sin(y) cos(y))
            /(cos^2(y) cosh^2(x) + sin^2(y) sinh^2(x))
--]]

function math.tanh(n)
    if type(n) == "number" then return __math.tanh(n) end
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            local x = n:real()
            local y = n:imaginary()
            local cy = __math.cos(y)
            local sy = __math.sin(y)
            local chx = __math.cosh(x)
            local shx = __math.sinh(x)
            local d = cy^2 * chx^2 + sy^2 * shx^2
            if d == 0 then
                return false
            end
            return Complex(shx*chx/d,sy*cy/d)
    end
    print("Cannot take the hyperbolic tangent of " .. n)
end

--[[
log(r e^(i a)) = log(r) + i (a + 2 k pi)
--]]

function math.log(n,k)
    if type(n) == "number" then
        if k then
            return Complex(__math.log(n), 2*k*math.pi)
        else
            return __math.log(n)
        end
    end
    k = k or 0
    if type(n) == "table"
        and n.is_a
        and n:is_a(Complex)
        then
            return Complex(__math.log(n:len()),n:arg() + 2*k*math.pi)
    end
    print("Cannot take the logarithm of " .. n)
end

if cmodule.loaded "TestSuite" then
    
    testsuite.addTest({
        name = "Complex",
        setup = function()
    local z = Complex(2,4)
    local w = Complex(-1,1)
    for k,v in ipairs({
        {"z", z},
        {"w", w},
        {"Muliplication" , z*w},
        {"ScalingLeft" , 2*z},
        {"ScalingRight" , z*2},
        {"Addition" , z+w},
        {"Addition" , 3+w},
        {"Addition" , w+3},
        {"Subtraction" , z-w},
        {"Subtraction" , 3-w},
        {"Subtraction" , w-3},
        {"Division" , z/w},
        {"DivisionLeft" , 2/z},
        {"DivisionRight" , z/2},
        {"Reciprocal", 1/z},
        {"Negation", -z},
        {"Powers", z^2},
        {"Roots", z^.5},
        {"Complex Powers", z^w},
        {"Conjugation", w^""},
        {"Polar Form", z:topolarstring()},
        {"Length", z:len()},
        {"Overloading",""},
        {"Absolute value", math.abs(z)},
        {"Powers", math.pow(z,w)},
        {"Square Root", math.sqrt(z)},
        {"Cosine", math.cos(z)},
        {"Sine", math.sin(z)},
        {"Tangent", math.tan(z)},
        {"Hyperbolic Cosine", math.cosh(z)},
        {"Hyperbolic Sine", math.sinh(z)},
        {"Hyperbolic Tangent", math.tanh(z)},
        {"Logarithm", math.log(z,1)},
    }) do
        print(v[1] .. ": " .. v[2])
    end
end,
draw = function()
end
})

end

return Complex

--]==]

Main.lua

--Project: Library Maths
--Version: 2.1
--Dependencies:
--Comments: 

VERSION = 2.1
clearProjectData()
-- DEBUG = true
-- Use this function to perform your initial setup
function setup()
    if AutoGist then
        autogist = AutoGist("Library Maths","A library of classes and functions for mathematical objects.",VERSION)
        autogist:backup(true)
    end
    if not cmodule then
        openURL("http://loopspace.mathforge.org/discussion/36/my-codea-libraries")
        print("You need to enable the module loading mechanism.")
        print("See http://loopspace.mathforge.org/discussion/36/my-codea-libraries")
        print("Touch the screen to exit the program.")
        draw = function()
        end
        touched = function()
            close()
        end
        return
    end
    --displayMode(FULLSCREEN_NO_BUTTONS)
    cmodule "Library Maths"
    cmodule.path("Library Base", "Library UI", "Library Utilities")
    cimport "TestSuite"
    Complex = cimport "Complex"
    Quaternion = cimport "Quaternion"
    --displayMode(FULLSCREEN)
    local Touches = cimport "Touch"
    local UI = cimport "UI"
    cimport "Menu"
    touches = Touches()
    ui = UI(touches)
    ui:systemmenu()
    fps = {}
    for k=1,20 do
        table.insert(fps,1/60)
    end
    afps = 60
    parameter.watch("math.floor(20/afps)")
    sl = Quaternion.unit():make_slerp(Quaternion(1.00001,0,0,0))
    print(sl(1))
    qc = vec3(1,0,0):rotateTo(vec3(0,0,1))
    print(qc)
    print(qc:tomatrix())
    print(vec3(1,0,0)^qc)
    print(vec3(0,1,0)^qc)
    print(vec3(0,0,1)^qc)
    print(qc*Quaternion(0,1,0,0)*qc^"")
    print(qc^""*Quaternion(0,1,0,0)*qc)
    qq = Quaternion(qc)
end

function draw()
    touches:draw()
    table.remove(fps,1)
    table.insert(fps,DeltaTime)
    afps = 0
    for k,v in ipairs(fps) do
        afps = afps + v
    end
    background(75, 104, 90, 255)
    ui:draw()
end

function touched(touch)
    touches:addTouch(touch)
end

function orientationChanged(o)
end

function fullscreen()
end

function reset()
end

Matrix.lua

--[==[
-- Matrix
local Vector = cimport "Vector"
local _,Sentence = unpack(cimport "Font")
local Matrix = class()

function Matrix:init(t)
    local r = 0
    local c = 0
    for k,v in ipairs(t) do
        table.insert(self,Vector(v))
        r = r + 1
    end
    c = self[1].size
    for k,v in ipairs(self) do
        if v.size ~= c then
            return false
        end
    end
    
    self.rows = r
    self.cols = c
end

function Matrix:add(m)
    if self.rows ~= m.rows then
        return false
    end
    if self.cols ~= m.cols then
        return false
    end
    local u = {}
    for k,v in ipairs(self) do
        table.insert(u,v + m[k])
    end
    return Matrix(u)
end

function Matrix:subtract(m)
    if self.rows ~= m.rows then
        return false
    end
    if self.cols ~= m.cols then
        return false
    end
    local u = {}
    for k,v in ipairs(self) do
        table.insert(u,v - m[k])
    end
    return Matrix(u)
end

function Matrix:scale(l)
    local u = {}
    for k,v in ipairs(self) do
        table.insert(u,l*v)
    end
    return Matrix(u)
end

function Matrix:multiplyRight(m)
    if self.cols ~= m.rows then
        return false
    end
    local u = {}
    for k,v in ipairs(self) do
        table.insert(u,v*m)
    end
    return Matrix(u)
end

function Matrix:multiplyLeft(m)
    return m:multiplyRight(self)
end

function Matrix:transpose()
    local nrows = {}
    for i = 1,self.cols do
        nrows[i] = {}
    end
    for k,v in ipairs(self) do
        for i = 1,self.cols do
            table.insert(nrows[i],v[i])
        end
    end
    return Matrix(nrows)
end

--[[
Inline operators:

u + v
u - v
-u
u * v : scaling
u / v : scaling
u == v : equality
--]]

function Matrix:__add(v)
    return self:add(v)
end

function Matrix:__sub(v)
    return self:subtract(v)
end

function Matrix:__unm()
    return self:scale(-1)
end

function Matrix:__mul(v)
    if type(self) == "number" then
        return v:scale(self)
    elseif type(v) == "number" then
        return self:scale(v)
    elseif type(self) == "userdata" then
        return Vector(self):applyMatrixRight(v):tovec()
    elseif type(v) == "userdata" then
        return Vector(v):applyMatrixLeft(self):tovec()
    else
        if self.is_a and self:is_a(Vector) then
            return self:applyMatrixRight(v)
        elseif v.is_a and v:is_a(Vector) then
            return v:applyMatrixLeft(self)
        else
            return self:multiplyRight(v)
        end
    end
    return false
end

function Matrix:__div(l)
    if type(l) == "number" then
        return self:scale(1/l)
    else
        return false
    end
end

function Matrix:__eq(v)
    return self:is_eq(v)
end

function Matrix:__concat(v)
    if type(v) == "table" 
        and v:is_a(Matrix) then
            return self .. v:tostring()
        else
            return self:tostring() .. v
        end
end

function Matrix:tostring()
    local u = {}
    for k,v in ipairs(self) do
        table.insert(u,v:tostring())
    end
    return "[" .. table.concat(u,";") .. "]"
end

function Matrix:__tostring()
    return self:tostring()
end

function Matrix:toModelMatrix()
    local m = {}
    for j=1,4 do
        for i=1,4 do
            if self[i] and self[i][j] then
                table.insert(m,self[i][j])
            elseif i == j then
                table.insert(m,1)
            else
                table.insert(m,0)
            end
        end
    end
    return matrix(unpack(m))
end

function Matrix:setParameters(t)
    t = t or {}
    for k,v in pairs({
        "font",
        "colour",
        "ui",
        "pos",
        "rowSep",
        "colSep",
        "anchor"
    }) do
        if t[v] then
            self[v] = t[v]
        end
    end
    if t.touchHandler then
        t.touchHandler:pushHandler(self)
    end
end

function Matrix:prepare()
    if not self.prepared then
    local f = self.font
    local t = {}
    local r
    local w = {}
    local s
    local col = self.colour or Colour.svg.White
    for k,v in ipairs(self) do
        r = {}
        for l,u in ipairs(v) do
            s = Sentence(f,math.floor(100*u+.5)/100)
            s:prepare()
            if not w[l] then
                w[l] = s.width
            else
                w[l] = math.max(w[l],s.width)
            end
            s:setColour(col)
            table.insert(r,s)
        end
        table.insert(t,r)
    end
    self.entries = t
    self.widths = w
    self.prepared = true
    end
end

function Matrix:draw()
    pushStyle()
    resetStyle()
    self:prepare()
    local rs = self.rowSep or 0
    local cs = self.colSep or 0
    local x,y = WIDTH/2,HEIGHT/2
    if self.pos then
        x,y = self.pos()
    end
    local a = self.anchor or "centre"
    local widths = self.widths
    local w = cs * self.cols
    for k,v in ipairs(widths) do
        w = w + v
    end
    local lh = self.font:lineheight() + rs
    local h = self.rows * lh
    x,y = RectAnchorAt(x,y,w,h,a)
    local yy = y + h + rs/2
    local xx = x + cs/2
    y = yy + rs/2 + self.font.descent
    for k,v in ipairs(self.entries) do
        x = xx
        y = y - lh
        for l,u in ipairs(v) do
            u:draw(x,y)
            x = x + widths[l] + cs
        end
    end
    xx = xx - cs/2
    yy = yy - rs/2
    smooth()
    strokeWidth(4)
    local col = self.colour or Colour.svg.White
    stroke(col)
    lineCapMode(PROJECT)
    line(xx-5,yy,xx,yy)
    line(xx-5,yy,xx-5,yy-h)
    line(xx-5,yy-h,xx,yy-h)
    line(xx+w+5,yy,xx+w,yy)
    line(xx+w+5,yy,xx+w+5,yy-h)
    line(xx+w+5,yy-h,xx+w,yy-h)
    popStyle()
end
            
function Matrix:isTouchedBy(touch)
    local rs = self.rowSep or 0
    local cs = self.colSep or 0
    local x,y = WIDTH/2,HEIGHT/2
    if self.pos then
        x,y = self.pos()
    end
    local a = self.anchor or "centre"
    local widths = self.widths
    local w = cs * self.cols
    for k,v in ipairs(widths) do
        w = w + v
    end
    local lh = self.font:lineheight() + rs
    local h = self.rows * lh
    x,y = RectAnchorAt(x,y,w,h,a)
    if touch.x < x then
        return false
    end
    if touch.x > x + w then
        return false
    end
    if touch.y < y then
        return false
    end
    if touch.y > y + h then
        return false
    end
    return true
end

function Matrix:processTouches(g)
    if g.type.ended and g.type.tap then
        local t = g.touchesArr[1]
        local rs = self.rowSep or 0
    local cs = self.colSep or 0
    local x,y = WIDTH/2,HEIGHT/2
    if self.pos then
        x,y = self.pos()
    end
    local a = self.anchor or "centre"
    local lh = self.font:lineheight() + rs
    local h = self.rows * lh
        local widths = self.widths
    local w = cs * self.cols
    for k,v in ipairs(widths) do
        w = w + v
    end
    x,y = RectAnchorAt(x,y,w,h,a)
        local r = math.floor((y + h - t.touch.y)/lh) + 1
        
        w = x
        local c
        for k,v in ipairs(widths) do
            w = w + v + cs
            if t.touch.x < w then
                c = k
                break
            end
        end
        self.ui:getNumberSpinner({
            action = function(v)
                    self[r][c] = v
                    self.prepared = false
                    return true
                end,
            value = self[r][c]
        })
    end
    if g.updated then
        g:noted()
    end
    if g.type.ended then
        g:reset()
    end
end

return Matrix

--]==]

Path.lua

--[==[
-- Path object

local Path = class()

local PATH_STEP = 1
local PATH_RES = .05*.05
local PATH_RENDER = 0
local PATH_CREATE = 1
local PATH_POINTS = 2

function Path:init(t)
    t = t or {}
    self.style = t.style or {}
    self.lastPoint = vec2(0,0)
    self.elements = t.elements or {}
    if t.touchHandler then
        t.touchHandler:pushHandler(self)
    end
    
    self.touchpt = {}
    self.m = mesh()
end

function Path:deletePoints()
    self.elements = {}
    self.lastPoint = vec2(0,0)
end

function Path:draw()
    pushStyle()
    resetStyle()
    self:applyStyle(PATH_RENDER)
    self.m:draw()
    if self.edit then
        self:applyStyle(PATH_POINTS)
        local r = self.ptRadius
        local lpt
        pushStyle()
        noStroke()
        for k,v in ipairs(self.elements) do
            ellipse(v[2].x,v[2].y,r)
            if v[1] == "curve" then
                popStyle()
                line(lpt.x,lpt.y,v[2].x,v[2].y)
                line(v[3].x,v[3].y,v[4].x,v[4].y)
                pushStyle()
                noStroke()
                ellipse(v[3].x,v[3].y,r)
                ellipse(v[4].x,v[4].y,r)
                lpt = v[4]
            else
                lpt = v[2]
            end
        end
        popStyle()
    end
    popStyle()
end

function Path:generate()
    self:applyStyle(PATH_CREATE)
    local ver = {}
    local s,h
    for k,v in ipairs(self.elements) do
        ver,s,h = self:getPoints(v,ver,s,h)
    end
    if h then
        ver,s = HobbyPoints(ver,s,h)
    end
    self.lastPoint = s
    --debug:log({name = "path", message = "got called"})
    print("got called")
    self:makeMesh(ver)
end

function Path:applyStyle(t)
    local s = self.style or {}
    if t == PATH_CREATE then
        self.linewidth = s.linewidth or 5
        self.blur = s.blur or 1
        self.smooth = s.smooth or false
        self.drawColour = s.drawColour or Colour.svg.Black
    elseif t == PATH_RENDER then
    elseif t == PATH_POINTS then
        self.ptRadius = s.pointRadius or 15
        local l = s.pointLineWidth or 3
        strokeWidth(l)
        local l = s.pointLineCap or SQUARE
        lineCapMode(l)
    end
end

function Path:addElement(e)
    table.insert(self.elements,e)
end

function Path:moveto(v)
    self:addElement({"move",v})
    self.lastPoint = v
end

function Path:lineto(v)
    self:addElement({"line",v})
    self.lastPoint = v
end

function Path:curveto(b,c,d)
    self:addElement({"curve",b,c,d})
    self.lastPoint = d
end

function Path:curvethrough(v)
    self:addElement({"hobby",v})
    self.lastPoint = v
end

function Path:getPoints(e,ver,s,h)
    local t = e[1]
    if h and t ~= "hobby" then
        ver,s = HobbyPoints(ver,s,h)
    end
    if t == "move" then
        table.insert(ver,{e[2]})
        return ver, e[2]
    elseif t == "line" then
        local n = e[2] - s
        n = n:rotate90()
        n = n/n:len()
        table.insert(ver,{s,n})
        table.insert(ver,{e[2],n})
        return ver,e[2]
    elseif e[1] == "curve" then
        return BezierPoints(ver,s,e[2],e[3],e[4])
    elseif e[1] == "hobby" then
        h = h or {}
        table.insert(h,e[2])
        return ver,s,h
    else
        return ver,s
    end
end

function Path:makeMesh(pts)
    local ver = {}
    local col = {}
    local l = self.linewidth/2
    local b = self.blur + l
    local s = false -- self.smooth
    local c = self.drawColour
    local ct = Colour.opacity(c,0)
    local p,n
    local m = 0
    for k,v in ipairs(pts) do
        if v[2] and n then
        table.insert(ver,p + l*n)
        table.insert(ver,p - l*n)
        table.insert(ver,v[1] + l*v[2])
        table.insert(ver,v[1] + l*v[2])
        table.insert(ver,v[1] - l*v[2])
        table.insert(ver,p - l*n)
        for i=1,6 do
            table.insert(col,c)
        end
    m = m + 6
        if s then
        table.insert(ver,p + l*n)
        table.insert(ver,p + b*n)
        table.insert(ver,v[1] + l*v[2])
        table.insert(ver,v[1] + l*v[2])
        table.insert(ver,v[1] + b*v[2])
        table.insert(ver,p + b*n)
        table.insert(ver,p - l*n)
        table.insert(ver,p - b*n)
        table.insert(ver,v[1] - l*v[2])
        table.insert(ver,v[1] - l*v[2])
        table.insert(ver,v[1] - b*v[2])
        table.insert(ver,p - b*n)
        for i=1,12 do
            table.insert(col,ct)
        end
        end
        end
        p,n = unpack(v)
        
    end
    self.m.vertices = ver
    self.m.colors = col
end

function Path:isTouchedBy(touch)
    local p = vec2(touch.x,touch.y)
    for k,v in ipairs(self.elements) do
        if p:distSqr(v[2]) < 625 then
            self.touchpt[touch.id] = v[2]
            self.edit = true
            return true
        end
        if v[1] == "curve" then
            if p:distSqr(v[3]) < 625 then
                self.touchpt[touch.id] = v[3]
                self.edit = true
                return true
            end
            if p:distSqr(v[4]) < 625 then
                self.touchpt[touch.id] = v[4]
                self.edit = true
                return true
            end
        end
    end
    self.edit = false
    return false
end

function Path:processTouches(g)
    local regenerate = false
    if g.updated then
        for k,t in ipairs(g.touchesArr) do
            if t.updated then
                regenerate = true
                self.touchpt[t.touch.id].x = 
                    self.touchpt[t.touch.id].x
                    + t.touch.deltaX
                self.touchpt[t.touch.id].y = 
                    self.touchpt[t.touch.id].y
                    + t.touch.deltaY
            end
            if t.touch.state == ENDED then
                self.touchpt[t.touch.id] = nil
            end
        end
        g:noted()
        if regenerate then
            self:generate()
        end
    end
    if g.type.ended then
        g:reset()
    end
end

function Path:setLineWidth(l)
    self.style.linewidth = l
end

function Path:setBlur(l)
    self.style.blur = l
end

function Path:setSmooth(l)
    self.style.smooth = l
end

function Path:setColour(c)
    self.style.drawColour = c
end

function BezierPoints(pts,a,b,c,d)
    local t = 0
    local r = PATH_RES
    local s = PATH_STEP
    local tpt = a
    table.insert(pts,{tpt,Bezierunml(0,a,b,c,d)})
        local dis
        local p
    while t < 1 do
        dis = 0
        while dis < r do
            t = t + s
            p = Bezierpt(t,a,b,c,d)
            dis = p:distSqr(tpt)
        end
        if t > 1 then
            t = 1
            p = d
        end
        table.insert(pts,{p,Bezierunml(t,a,b,c,d)})
        tpt = p
    end
    return pts,d
    
end

function Bezierpt(t,a,b,c,d)
    local s = 1 - t
    return s^3 * a
        + 3*s*s*t * b
        + 3*s*t*t * c
        + t^3 * d
end

function Beziertgt(t,a,b,c,d)
    local s = 1 - t
    return 3*s^2 * (b - a)
        + 6*s*t * (c - b)
        + 3*t^2 * (d - c)
end

function Beziernml(t,a,b,c,d)
    return Beziertgt(t,a,b,c,d):rotate90()
end

function Bezierunml(t,a,b,c,d)
    local pt = Beziernml(t,a,b,c,d)
    local l = pt:len()
    if l == 0 then
        return vec2(0,0)
    else
        return pt/l
    end
end

function Bezierutgt(t,a,b,c,d)
    local pt = Beziertgt(t,a,b,c,d)
    local l = pt:len()
    if l == 0 then
        return vec2(0,0)
    else
        return pt/l
    end
end

function QuickHobby(a,b,c,tha)
    if type(a) == "table" then
        tha = b
        a,b,c = unpack(a)
    end
    local da = a:dist(b)
    local db = b:dist(c)
    local wa = vec2(1,0):angleBetween(b-a)
    local wb = vec2(1,0):angleBetween(c-b)
    local psi = wb - wa
    if psi > math.pi then
        psi = psi - 2*math.pi
    end
    if psi <= -math.pi then
        psi = psi + 2*math.pi
    end
    local thb,phb,phc
    if tha then
        thb = -(2*psi + tha) * db / (2*db + da)
        phb = - psi - thb
        phc = thb
    else
        thb = - psi * db / (da + db)
        tha = - psi - thb
        phb = tha
        phc = thb
    end
    local alpha = math.sqrt(2) * (math.sin(tha) - 1/16 * math.sin(phb)) * (math.sin(phb) - 1/16 * math.sin(tha)) * (math.cos(tha) - math.cos(phb))
    local rho = (2 + alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(tha) + (3 - math.sqrt(5))/2 * math.cos(phb))
    local sigma = (2 - alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(phb) + (3 - math.sqrt(5))/2 * math.cos(tha))
    local ctrlaa = a + da * rho * vec2(math.cos(tha + wa), math.sin(tha + wa))/3
    local ctrlab = b - da * sigma * vec2(math.cos(wa - phb), math.sin(wa - phb))/3
    alpha = math.sqrt(2) * (math.sin(thb) - 1/16 * math.sin(phc)) * (math.sin(phc) - 1/16 * math.sin(thb)) * (math.cos(thb) - math.cos(phc))
    rho = (2 + alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(thb) + (3 - math.sqrt(5))/2 * math.cos(phc))
    sigma = (2 - alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(phc) + (3 - math.sqrt(5))/2 * math.cos(thb))
    local ctrlba = b + db * rho * vec2(math.cos(thb + wb), math.sin(thb + wb))/3
    local ctrlbb = c - db * sigma * vec2(math.cos(wb - phc), math.sin(wb - phc))/3
    return {a,ctrlaa,ctrlab,b},{b,ctrlba,ctrlbb,c},thb
end

function QHobbyGenerator(a,b)
    local th
    return function(c)
        local p,q
        p,q,th = QuickHobby(a,b,c,th)
        a = b
        b = c
        return p,q
    end
end

function HobbyPoints(ver,s,pts,extra)
    if #pts == 1 then
        local v = pts[1]
        if type(v) == "table" then
            v = v[1]
        end
        local n = v - s
        n = n:rotate90()
        n = n/n:len()
        table.insert(ver,{s,n})
        table.insert(ver,{v,n})
        return ver,v
    end
    table.insert(pts,1,s)
    local bcurves = Hobby(pts,extra)
    for _,v in ipairs(bcurves) do
        ver,s = BezierPoints(ver,unpack(v))
    end
    return ver,s
end

function Hobby(pts,extra)
    local z = {}
    local d = {}
    local omega = {}
    local psi = {}
    local it = {}
    local ot = {}
    local n = 0
    local A = {}
    local B = {}
    local C = {}
    local D = {}
    local icurl = 1
    local ocurl = 1
    local theta = {}
    local phi = {}
    local rho ={}
    local sigma = {}
    local a = math.sqrt(2)
    local b = 1/16
    local c = (3 - math.sqrt(5))/2
    local cpta = {}
    local cptb = {}
    local dten = 1
    local curves = {}
    ot[0] = 1
    if extra then
        dten = extra.tension or 1
        icurl = extra.inCurl or 1
        ocurl = extra.outCurl or 1
        ot[0] = extra.initialTension or 1
    end

    local ang,i,o,v
    z[0] = table.remove(pts,1)
    for _,t in ipairs(pts) do
        if type(t) == "table" then
            v = t[1]
            i = t[2] or dten
            o = t[3] or dten
        else
            v = t
            i = dten
            o = dten
        end
        d[n] = z[n]:dist(v)
        omega[n] = vec2(1,0):angleBetween(v-z[n])
       if n > 0 then
          ang = omega[n] - omega[n-1]
       if ang > math.pi then
             ang = ang - 2*math.pi
        end
        if ang < -math.pi then
        ang = ang + 2*math.pi
         end
          psi[n] = ang
       end
       n = n + 1
       z[n] = v
       it[n] = i
       ot[n] = o
    end

    if extra then
        theta[0] = extra.inAngle
        phi[n] = extra.outAngle
    end

   for k=1,n-1 do
      A[k] = d[k] * it[k+1] * it[k]^2
   end

   if theta[0] then
      B[0] = 1
   else
      B[0] = ot[0]^3 * (3 * it[1] - 1) + icurl * it[1]^3
   end
   for k = 1,n-2 do
      B[k] = d[k] * it[k+1] * it[k]^2 * (3 * ot[k-1] - 1) + d[k-1] * ot[k-1] * ot[k]^2 * (3 * it[k+1] - 1)
   end
   B[n-1] = d[n-1] * it[n] * it[n-1]^2 * (3 * ot[n-2] - 1) + d[n-2] * ot[n-2] * ot[n-1]^2 * (3 * it[n] - 1)
   if not phi[n] then
      B[n-1] = B[n-1] - d[n-2] * ot[n-2] * ot[n-1]^2 * (it[n]^3 + ocurl * ot[n-1]^3 * (3 * it[n] - 1)) / (it[n]^3 * (3 * ot[n-1] - 1) + ocurl * ot[n-1]^3)
   end

   if theta[0] then
      C[0] = 0
   else
      C[0] = ot[0]^3 + icurl * it[1]^3 * (3 * ot[0] - 1)
   end
   for i=1,n do
      C[i] = d[i-1] * ot[i-1] * ot[i]^2
   end

   if theta[0] then
      D[0] = theta[0]
   else
      D[0] = - (ot[0]^3 + icurl * it[1]^3 * (3 * ot[0] - 1)) * psi[1]
   end
   for i=1,n-2 do
      D[i] = - d[i] * it[i+1] * it[i]^2 * (3 * ot[i-1] - 1) * psi[i] - d[i-1] * ot[i-1] * ot[i]^2 * psi[i+1]
   end
   D[n-1] = - d[n-1] * it[n] * it[n-1]^2 * (3 * ot[n-2] - 1) * psi[n-1]
   if phi[n] then
      D[n-1] = D[n-1] - d[n-2] * ot[n-2] * ot[n-1]^2 * phi[n]
   end

   for i=1,n-1 do
      B[i] = B[i-1] * B[i] - A[i] * C[i-1]
      C[i] = B[i-1] * C[i]
      D[i] = B[i-1] * D[i] - A[i] * D[i-1]
   end
   theta[n-1] = D[n-1]/B[n-1]
   for i=n-2,1,-1 do
      theta[i] = (D[i] - C[i] * theta[i+1])/B[i]
   end
   for i=1,n-1 do
      phi[i] = -psi[i] - theta[i]
   end
   if not theta[0] then
      theta[0] = (ot[0]^3 + icurl * it[1]^3 * (3 * ot[0] - 1)) / (ot[0]^3 * (3 * it[1] - 1) + icurl * it[1]^3) * phi[1]
   end
   if not phi[n] then
      phi[n] = (it[n]^3 + ocurl * it[n-1]^3 * (3 * it[n] - 1)) / (it[n]^3 * (3 * ot[n-1] - 1) + ocurl * ot[n-1]^3) * theta[n-1]
   end

   local alpha
   for i = 0,n-1 do
      alpha = a * (math.sin(theta[i]) - b * math.sin(phi[i+1])) * (math.sin(phi[i+1]) - b * math.sin(theta[i])) * (math.cos(theta[i]) - math.cos(phi[i+1]))
      rho[i] = (2 + alpha) / (1 + (1 - c) * math.cos(theta[i]) + c * math.cos(phi[i+1]))
      sigma[i+1] = (2 - alpha) / (1 + (1 - c) * math.cos(phi[i+1]) + c * math.cos(theta[i]))
   end

   for i = 0,n-1 do
    table.insert(curves,{
        z[i],
        z[i] + d[i]*rho[i] * vec2(math.cos(theta[i] + omega[i]), math.sin(theta[i] + omega[i]))/3,
        z[i+1] - d[i] * sigma[i+1] * vec2(math.cos(omega[i] - phi[i+1]), math.sin(omega[i] - phi[i+1]))/3,
        z[i+1]
        })
   end
    return curves
end

cmodule.gexport {
   BezierPoints = BezierPoints,
   Bezierpt = Bezierpt,
   Beziertgt = Beziertgt,
   Beziernml = Beziernml,
   Bezierunml = Bezierunml,
   Bezierutgt = Bezierutgt,
   QuickHobby = QuickHobby,
   HobbyPoints = HobbyPoints,
    QHobbyGenerator = QHobbyGenerator,
    Hobby = Hobby
	}

return Path

--]==]

Quaternion.lua

--[==[
-- 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 Quaternion = class()

--[[
A quaternion can either be specified by giving the four coordinates as
real numbers or by giving the scalar part and the vector part.
--]]

local error = error or print

function Quaternion:init(...)
    -- you can accept and set parameters here
    if arg.n == 4 then
        -- four numbers
        self.q = vec4(arg[1],arg[2],arg[3],arg[4])
    elseif arg.n == 2 then
        -- real number plus vector
        self.q = vec4(arg[1],arg[2].x,arg[2].y,arg[2].z)
    elseif arg.n == 1 then
        if type(arg[1]) == "userdata" then
            self.q = vec4(arg[1].x,arg[1].y,arg[1].z,arg[1].w)
        elseif type(arg[1]) == "table" and arg[1].is_a and arg[1]:is_a(Quaternion) then
            self.q = vec4(arg[1].q.x,arg[1].q.y,arg[1].q.z,arg[1].q.w)
        else
            error("Incorrect type of argument to Quaternion")
        end
    else
        error("Incorrect number of arguments to Quaternion")
    end
end

--[[
Test if we are zero.
--]]

function Quaternion:is_zero()
    return self.q == vec4(0,0,0,0)
end

--[[
Test if we are real.
--]]

function Quaternion:is_real()
    -- are we the zero vector
    if self.q.y ~= 0 or self.q.z ~= 0 or self.q.w ~= 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.q.x ~= 0 then
        return false
    end
    return true
end

--[[
Test for equality.
--]]

function Quaternion:is_eq(q)
    return self.q == q.q
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.q:dot(q.q)
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 self.q:len()
end

--[[
Often enough to know the length squared, which is quicker.
--]]

function Quaternion:lenSqr()
    return self.q:lenSqr()
end

--[[
Distance between two quaternions.
--]]

function Quaternion:dist(q)
    return self.q:dist(q.q)
end

--[[
Often enough to know the distance squared, which is quicker.
--]]

function Quaternion:distSqr(q)
    return self.q:distSqr(q.q)
end

--[[
Normalise a quaternion to have length 1, if possible.
--]]

function Quaternion:normalise()
    return Quaternion(self.q:normalize())
end

function Quaternion:normalize()
    return Quaternion(self.q:normalize())
end

--[[
Scale the quaternion.
--]]

function Quaternion:scale(l)
    return Quaternion(self.q * l)
end

--[[
Add two quaternions.
--]]

function Quaternion:add(q)
    if type(q) == "number" then
        return Quaternion(self.q.x + q, self.q.y, self.q.z, self.q.w)
    else
        return Quaternion(self.q + q.q)
    end
end

--[[
q + p
--]]

function Quaternion:__add(q)
    if type(q) == "number" then
        return Quaternion(self.q.x + q, self.q.y, self.q.z, self.q.w)
    elseif type(self) == "number" then
        return Quaternion(self + q.q.x, q.q.y, q.q.z, q.q.w)
    else
        return Quaternion(self.q + q.q)
    end
end

--[[
Subtraction
--]]

function Quaternion:subtract(q)
    return Quaternion(self.q - q.q)
end

--[[
q - p
--]]

function Quaternion:__sub(q)
    if type(q) == "number" then
        return Quaternion(self.q.x - q, self.q.y, self.q.z, self.q.w)
    elseif type(self) == "number" then
        return Quaternion(self - q.q.x, - q.q.y, - q.q.z, - q.q.w)
    else
        return Quaternion(self.q - q.q)
    end
end

--[[
Negation (-q)
--]]

function Quaternion:__unm()
    return Quaternion(-self.q)
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.q.x * q.q.x - self.q.y * q.q.y - self.q.z * q.q.z - self.q.w * q.q.w
    b = self.q.x * q.q.y + self.q.y * q.q.x + self.q.z * q.q.w - self.q.w * q.q.z
    c = self.q.x * q.q.z - self.q.y * q.q.w + self.q.z * q.q.x + self.q.w * q.q.y
    d = self.q.x * q.q.w + self.q.y * q.q.z - self.q.z * q.q.y + self.q.w * q.q.x
    return Quaternion(a,b,c,d)
end

--[[
q * p
--]]

function Quaternion:__mul(q)
    if type(q) == "number" then
        return self:scale(q)
    elseif type(self) == "number" then
        return q:scale(self)
    elseif type(q) == "table" and q.is_a and q:is_a(Quaternion) 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.q.x, - self.q.y, - self.q.z, - self.q.w)
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:lenSqr()
    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)
    elseif type(n) == "table" and n.is_a and n:is_a(Quaternion) then
        return self:multiplyLeft(n):multiplyRight(n:reciprocal())
    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(self) == "number" then
        return q:reciprocal():scale(self)
    elseif type(q) == "table" and q.is_a and q:is_a(Quaternion) then
        return self:multiplyRight(q:reciprocal())
    end
end


--[[
Interpolation functions, we assume the input to be already normalised
for speed.  If you cannot guarantee this, renormalise the input first.
The constructor functions do do the renormalisation.
--]]

--[[
Linear interpolation, renormalised
--]]

function Quaternion:lerp(q,t)
    if not t then
        return Quaternion.unit():lerp(self,q)
    end
    local v
    if self.q == -q.q then
        -- antipodal points, need a midpoint
        v = vec4(self.q.y,-self.q.x,self.q.w,-self.q.z)
        v = (1 - 2*t)*self.q + (1-math.abs(2*t-1))*v
    else
        v = (1-t)*self.q + t*q.q
    end
    return Quaternion(v:normalize())
end

--[[
Spherical interpolation
--]]

function Quaternion:slerp(q,t)
    if not t then
        return Quaternion.unit():slerp(self,q)
    end
    local v
    if self.q == -q.q then
        -- antipodal points, need a midpoint
        v = vec4(self.q.y,-self.q.x,self.q.w,-self.q.z)
        t = 2*t
    elseif self.q == q.q then
        return Quaternion(self)
    else
        v = q.q
    end
    local ca = self.q:dot(v)
    local sa = math.sqrt(1 - math.pow(ca,2))
    if sa == 0 then
        return Quaternion(self)
    end
    local a = math.acos(ca)
    sa = math.sin(a*t)/sa
    v = (math.cos(a*t) - ca*sa)*self.q+ sa*v
    return Quaternion(v)
end

--[[
Constructor for normalised linear interpolation.
--]]

function Quaternion:make_lerp(q)
    if not q then
        return Quaternion.unit():make_lerp(self)
    end
    local v,w
    w = self.q:normalize()
    if self.q == -q.q then
        -- antipodal points, need a midpoint
        v = vec4(w.y,-w.x,w.w,-w.z)
        return function(t)
            local u = (1 - 2*t)*w + (1-math.abs(2*t-1))*v
            return Quaternion(u:normalize())
        end
    else
        v = q.q:normalize()
        return function(t)
            local u = (1-t)*w + t*v
            return Quaternion(u:normalize())
        end
    end
end

--[[
Spherical interpolation
--]]

function Quaternion:make_slerp(q)
    if not q then
        return Quaternion.unit():make_slerp(self)
    end
    local v,f,u
    if self.q == -q.q then
        -- antipodal points, need a midpoint
        v = vec4(self.q.y,-self.q.x,self.q.w,-self.q.z)
        f = 2
    elseif self.q == q.q then
        return function(t)
            return Quaternion(self)
        end
    else
        v = q.q
        f = 1
    end
    v = v:normalize()
    u = self.q:normalize()
    local ca = u:dot(v)
    local sa = math.sqrt(1 - math.pow(ca,2))
    if sa == 0 then
        return function(t)
            return Quaternion(self)
        end
    end
    local a = math.acos(ca)
    v = (v - ca*self.q)/sa
    return function(t)
        local u = math.cos(a*f*t)*self.q + math.sin(a*f*t)*v
        return Quaternion(u)
    end
end

--[[
Returns the real part.
--]]

function Quaternion:real()
    return self.q.x
end

--[[
Returns the vector (imaginary) part as a Vec3 object.
--]]

function Quaternion:vector()
    return vec3(self.q.y, self.q.z, self.q.w)
end

--[[
Represents a quaternion as a string.
--]]

function Quaternion:tostring()
    local s
    local im ={{self.q.y,"i"},{self.q.z,"j"},{self.q.w,"k"}}
    if self.q.x ~= 0 then
        s = string.format("%.3f",self.q.x)
    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 .. " + " .. 
                        string.format("%.3f",v[1]) .. v[2]
                    end
                else
                    if v[1] == -1 then
                        s = s .. " - " .. v[2]
                    else
                        s = s .. " - " .. 
                        string.format("%.3f",-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 = string.format("%.3f",v[1]) .. v[2]
                end
        end
    end
    end
    if s then
        return s
    else
        return "0"
    end
end

function Quaternion:__tostring()
    return self:tostring()
end

function Quaternion:__concat(s)
    if type(s) == "string" then
        return self:tostring() .. s
    elseif type(s) == "table" and s.is_a and s:is_a(Quaternion) then
        return self .. s:tostring()
    end
end

function Quaternion:tomatrix()
    local a,b,c,d = self.q.x,-self.q.y,-self.q.z,-self.q.w
    local ab = 2*a*b
    local ac = 2*a*c
    local ad = 2*a*d
    local bb = 2*b*b
    local bc = 2*b*c
    local bd = 2*b*d
    local cc = 2*c*c
    local cd = 2*c*d
    local dd = 2*d*d
    return matrix(
        1 - cc - dd,
        bc - ad,
        ac + bd,
        0,
        bc + ad,
        1 - bb - dd,
        cd - ab,
        0,
        bd - ac,
        cd + ab,
        1 - bb - cc,
        0,
        0,0,0,1
    )
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
--]]

function Quaternion:Gravity()
    local qg,qx
    if not self then
        qx = Quaternion.ReferenceFrame()
    else
        qx = self
    end
    local y = vec3(0,-1,0)^qx
    qg = y:rotateTo(Gravity)
    return qg*qx
end

local frame = {
        vec3(1,0,0),
        vec3(0,1,0),
        vec3(0,0,1)
    }
local qzyx = Quaternion(1,0,0,0)

--[[
Needs to be run once every frame!
--]]

function Quaternion:updateReferenceFrame()
    local x,y,z
    if CurrentOrientation == PORTRAIT then
        x,y,z = -frame[1],-frame[2],-frame[3]
    elseif CurrentOrientation == PORTRAIT_UPSIDE_DOWN then
        x,y,z = frame[1],frame[2],-frame[3]
    elseif CurrentOrientation == LANDSCAPE_LEFT then
        x,y,z = frame[2],-frame[1],-frame[3]
    elseif CurrentOrientation == LANDSCAPE_RIGHT then
        x,y,z = -frame[2],frame[1],-frame[3]
    end
    local qz = Quaternion.Rotation(RotationRate.z*DeltaTime,z.x,z.y,z.z)
    local qy = Quaternion.Rotation(RotationRate.y*DeltaTime,y.x,y.y,y.z)
    local qx = Quaternion.Rotation(RotationRate.x*DeltaTime,x.x,x.y,x.z)
    if self then
        local q = qz * qy * qx * self
        self.q = q.q
    else
        qzyx = qz * qy * qx * qzyx
        return qzyx
    end
end

function Quaternion:ReferenceFrame()
    return self or qzyx
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:normalise()
    c = math.cos(a/2)
    s = math.sin(a/2)
    q = q:scale(s)
    q = q:add(c)
    return q
end

--[[
The unit quaternion.
--]]

function Quaternion.unit()
    return Quaternion(1,0,0,0)
end

--[[
Extensions to vec3 type.
--]]

do
    local mt = getmetatable(vec3())

--[[
Promote to a quaternion with 0 real part.
--]]

mt["toQuaternion"] = function (self)
    return Quaternion(0,self.x,self.y,self.z)
end

--[[
Apply a quaternion as a rotation (assumes unit quaternion for speed).
--]]

mt["applyQuaternion"] = function (self,q)
    local x = self:toQuaternion()
    x = q:multiplyRight(x)
    x = x:multiplyRight(q:conjugate())
    return x:vector()
end

mt["__pow"] = function (self,q)
    if type(q) == "table" then
        if q:is_a(Quaternion) then
                return self:applyQuaternion(q)
        end
    end
    return false
end

mt["rotateTo"] = function (self,v)
    if v:lenSqr() == 0 or self:lenSqr() == 0 then
        return Quaternion(1,0,0,0)
    end
    local u = self:normalize()
    v = u + v:normalize()
    if v:lenSqr() == 0 then
        -- Opposite vectors, no canonical direction
        local a,b,c = math.abs(u.x), math.abs(u.y), math.abs(u.z)
        if a < b and a < c then
            v = vec3(0,-u.z,u.y)
        elseif b < c then
            v = vec3(u.z,0,-u.x)
        else
            v = vec3(u.y,-u.x,0)
        end
    end
    v = v:normalize()
    return Quaternion(u:dot(v),u:cross(v))
end

end

if cmodule.loaded "TestSuite" then
    
    testsuite.addTest({
        name = "Quaternion",
        setup = function()
            local q = Quaternion(.3,.4,.5,.6)
            local qq = Quaternion(.2,vec3(.4,-.5,.1))
            for k,v in ipairs({
                {"Quaternion", q},
                {"Sum", q + qq},
                {"Sum (number)", q + 5},
                {"Sum (number)", 5 + q},
                {"Subtract", q - qq},
                {"Subtract (number)", q - 5},
                {"Subtract (number)", 5 - q},
                {"Product", q * qq},
                {"Scale (right)", q * 5},
                {"Scale (left)", 5 * q},
                {"Divison", q / qq},
                {"by number", q / 5},
                {"of number", 5 / q},
                {"Length", q:len()},
                {"Length Squared", q:lenSqr()},
            }) do
                print(v[1] .. ": " .. v[2])
            end
        end,
        draw = function()
        end
    })
end

return Quaternion

--]==]

Vec3.lua

--[==[
-- Vec3 class
-- Author: Andrew Stacey
-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
-- Licence: CC0 http://wiki.creativecommons.org/CC0

--[[
The "Vec3" class is for handling 3 dimensional vectors and defining a
variety of methods on them.

This is largely superceded by the vec3 userdata.  There are a few extra
methods, but these should really be added to vec3.
--]]

local Vec3 = class()

--[[
A 3-vector is three numbers.
--]]

function Vec3:init(x,y,z)
    self.x = x
    self.y = y
    self.z = z
end

--[[
Test for zero vector.
--]]

function Vec3:is_zero()
    if self.x ~= 0 or self.y ~= 0 or self.z ~= 0 then
        return false
    end
    return true
end

--[[
Test for equality.
--]]

function Vec3:is_eq(v)
    if self.x ~= v.x or self.y ~= v.y or self.z ~= v.z then
        return false
    end
    return true
end        

--[[
Inner product.
--]]

function Vec3:dot(v)
    return self.x * v.x + self.y * v.y + self.z * v.z
end

--[[
Cross product.
--]]

function Vec3:cross(v)
    local x,y,z
    x = self.y * v.z - self.z * v.y
    y = self.z * v.x - self.x * v.z
    z = self.x * v.y - self.y * v.x
    return Vec3(x,y,z)
end

--[[
Apply a given matrix (which is specified as a triple of vectors).
--]]

function Vec3:applyMatrix(a,b,c)
    local u,v,w
    u = a:scale(self.x)
    v = b:scale(self.y)
    w = c:scale(self.z)
    u = u:add(v)
    u = u:add(w)
    return u
end

--[[
Length of the vector
--]]

function Vec3:len()
    return math.sqrt(math.pow(self.x,2) + math.pow(self.y,2) + math.pow(self.z,2))
end

--[[
Squared length of the vector.
--]]

function Vec3:lenSqr()
    return math.pow(self.x,2) + math.pow(self.y,2) + math.pow(self.z,2)
end

--[[
Normalise the vector (if possible) to length 1.
--]]

function Vec3:normalise()
    local l
    if self:is_zero() then
        print("Unable to normalise a zero-length vector")
        return false
    end
    l = 1/self:len()
    return self:scale(l)
end

--[[
Scale the vector.
--]]

function Vec3:scale(l)
    return Vec3(self.x * l,self.y * l,self.z * l)
end

--[[
Add vectors.
--]]

function Vec3:add(v)
    return Vec3(self.x + v.x, self.y + v.y, self.z + v.z)
end

--[[
Subtract vectors.
--]]

function Vec3:subtract(v)
    return Vec3(self.x - v.x, self.y - v.y, self.z - v.z)
end

--[[
Apply a transformation between "absolute" coordinates and "relative"
coordinates.  In the "relative" system, xy are in the iPad screen and
z points straight out.  In the "absolute" system, y is in the
direction of the gravity vector, x is in the plane of the iPad screen,
and z is orthogonal to those two.

This function interprets the vector as being with respect to the
absolute coordinates and returns the corresponding vector in the
relative system.
--]]

function Vec3:absCoords()
    local gxy,l,ga,gb,gc
    gxy = Vec3(-Gravity.y,Gravity.x,0)
    l = gxy:len()
    if l == 0 then
    print("Unable to compute coordinate system, gravity vector is (" .. Gravity.x .. "," .. Gravity.y .. "," .. Gravity.z .. ")")
        return false
    end
    ga = gxy:scale(1/l)
    gb = Vec3(-Gravity.x,-Gravity.y,-Gravity.z)
    gc = Vec3(-Gravity.x * Gravity.z /l, -Gravity.y * Gravity.z /l, l)
    return self:applyMatrix(ga,gb,gc)
end

--[[
Determine whether or not the vector is in front of the "eye" (crude
test for visibility).
--]]

function Vec3:isInFront(e)
    if not e then
        e = Vec3.eye
    end
    if self:dot(e) < e:dot(e) then
        return true
    else
        return false
    end
end

--[[
Project the vector onto the screen using stereographic projection from
the "eye".
--]]

function Vec3:stereoProject(e)
    local t,v
    if not e then
        e = Vec3.eye
    end
    if self.z == e.z then
        -- can't project
        return false
    end
    t = 1 / (1 - self.z / e.z)
    v = self:subtract(e)
    v = v:scale(t)
    v = v:add(e)
    -- hopefully v.z is now 0!
    return vec2(v.x,v.y)
end

--[[
Partial inverse to stereographic projection: given a point on the
screen and a height, we find the point in space at that height which
would project to the point on the screen.
--]]

function Vec3.stereoInvProject(v,e,h)
    local t,u
    if not e then
        e = Vec3.eye
    end
    u = Vec3(v.x,v.y,0)
    t = h / e.z
    u = (1 - t) * u + t * e
    -- hopefully u.z is now h!
    return u
end

--[[
Returns the distance from the eye; useful for sorting objects.
--]]

function Vec3:stereoLevel(e)
    local v
    if not e then
        e = Vec3.eye
    end
    v = self:subtract(e)
    return e:len() - v:len()
end

--[[
Applies a rotation as specified by another 3-vector, with direction
being the axis and magnitude the angle.
--]]

function Vec3:rotate(w)
    local theta, u, v, a, b, c, x
    if w:is_zero() then
        return self
    else
        theta = w:len()
        w = w:normalise()
        if w.x ~= 0 then
                u = Vec3(-w.y/w.x,1,0)
                u = u:normalise()
        else
                u = Vec3(1,0,0)
        end
        v = w:cross(u)
        a = self:dot(u)
        b = self:dot(v)
        c = self:dot(w)
        x = w:scale(c)
        u = u:scale(a * math.cos(theta) + b * math.sin(theta))
        v = v:scale(-a * math.sin(theta) + b * math.cos(theta))
        x = x:add(u)
        x = x:add(v)
        return x
    end
end

--[[
Promote to a quaternion with 0 real part.
--]]

function Vec3:toQuaternion()
    return Quaternion(0,self.x,self.y,self.z)
end

--[[
Apply a quaternion as a rotation.
--]]

function Vec3:applyQuaternion(q)
    local x = self:toQuaternion()
    x = q:multiplyRight(x)
    x = x:multiplyRight(q:conjugate())
    return x:vector()
end

--[[
Inline operators:

u + v
u - v
-u
u * v : cross product (possibly bad choice) or scaling
u / v : scaling
u ^ q : apply quaternion as rotation
u == v : equality
u .. v : dot product (possibly bad choice)

The notation for cross product and dot product may be removed in a
later version.
--]]

function Vec3:__add(v)
    return self:add(v)
end

function Vec3:__sub(v)
    return self:subtract(v)
end

function Vec3:__unm()
    return self:scale(-1)
end

function Vec3:__mul(v)
    if type(self) == "number" then
        return v:scale(self)
    elseif type(v) == "number" then
        return self:scale(v)
    elseif type(v) == "table" then
        if v:is_a(Vec3) then
                return self:cross(v)
        end
    end
    return false
end

function Vec3:__div(l)
    if type(l) == "number" then
        return self:scale(1/l)
    else
        return false
    end
end

function Vec3:__pow(q)
    if type(q) == "table" then
        if q:is_a(Quaternion) then
                return self:applyQuaternion(q)
        end
    end
    return false
end

function Vec3:__eq(v)
    return self:is_eq(v)
end

function Vec3:__concat(v)
    if type(v) == "table" 
        and v:is_a(Vec3) 
        and type(self) == "table"
        and self:is_a(Vec3)
        then
        return self:dot(v)
    else
        if type(v) == "table" 
        and v:is_a(Vec3) then
            return self .. v:tostring()
        else
            return self:tostring() .. v
        end
    end
end

function Vec3:tostring()
    return "(" .. self.x .. "," .. self.y .. "," .. self.z .. ")"
end

function Vec3:__tostring()
    return self:tostring()
end

function Vec3:tovec3()
    return vec3(self.x,self.y,self.z)
end

--[[
The following functions are not class methods but are still related to
vectors.
--]]

--[[
Sets the "eye" for stereographic projection.  Input can either be a
Vec3 object or the information required to specify one.
--]]

function Vec3.SetEye(...)
    if arg.n == 1 then
        Vec3.eye = arg[1]
    elseif arg.n == 3 then
        Vec3.eye = Vec3(unpack(arg))
    else
        print("Wrong number of arguments to Vec3.SetEye (1 or 3 expected, got " .. arg.n .. ")")
        return false
    end
    return true
end

--[[
Some useful Vec3 objects.
--]]

Vec3.eye = Vec3(0,0,1)

Vec3.origin = Vec3(0,0,0)
Vec3.e1 = Vec3(1,0,0)
Vec3.e2 = Vec3(0,1,0)
Vec3.e3 = Vec3(0,0,1)

--[[
Is the line segment a-b over c-d when seen from e?
--]]

function Vec3.isOverLine(a,b,c,d,e)
    if not e then
        e = Vec3.eye
    end
    -- rebase at a
    b = b - a
    c = c - a
    d = d - a
    e = e - a
    -- test signs of various determinants
    a = c:cross(d)
    if a:dot(b) * a:dot(e) < 0 then
        return false
    end
    a = d:cross(e)
    if a:dot(b) * a:dot(c) < 0 then
        return false
    end
    a = e:cross(c)
    if a:dot(b) * a:dot(d) < 0 then
        return false
    end
    -- right direction, is it far enough?
    c = c:subtract(e)
    d = d:subtract(e)
    a = c:cross(d)
    local l = a:dot(b)
    local m = a:dot(e)
    if l * m > 0 and math.abs(l) > math.abs(m) then
        return true
    else
        return false
    end
end

--[[
Is the line segment a-b over the point c when seen from e?
r is the "significant distance"
--]]

function Vec3.isOverPoint(a,b,c,r,e)
    
    if not e then
        e = Vec3.eye
    end
    local aa,bb,ab,ac,bc,d,l
    -- rebase at e
    a = a - e
    b = b - e
    c = c - e
    d = a:cross(b)
    if math.abs(d:dot(c)) > r * d:len() then
        return false
    end
    aa = a:lenSqr()
    bb = b:lenSqr()
    ab = a:dot(b)
    ac = a:dot(c)
    bc = b:dot(c)
    
    if aa * bc < ab * ac then
        return false
    end
    if bb * ac < ab * bc then
        return false
    end
    
    l = math.sqrt((aa * bb - ab * ab) * (aa + bb - 2 * ab))
    
    if (bb - ab) * ac + (aa - ab) * bc < aa * bb - ab * ab + r * l then
        return false
    end
    return true
end

return Vec3

--]==]

Vector.lua

--[==[
-- Vector class
-- Author: Andrew Stacey
-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
-- Licence: CC0 http://wiki.creativecommons.org/CC0

--[[
The "Vector" class is for handling arbitrary dimensional vectors 
and defining a variety of methods on them.
--]]

local Vector = class()

--[[
A vector is a size and an array of numbers.
--]]

function Vector:init(...)
    local v
    if arg.n == 1 then
        if type(arg[1]) == "table" then
            if arg[1].is_a and arg[1]:is_a(Complex) then
                v = {arg[1]}
            else
                v = arg[1]
            end
        elseif type(arg[1]) == "userdata" then
            local mt = getmetatable(arg[1])
            if mt == getmetatable(vec2()) then
                v = {arg[1].x,arg[1].y}
            elseif mt == getmetatable(vec3()) then
                v = {arg[1].x,arg[1].y,arg[1].z}
            elseif mt == getmetatable(vec4()) then
                v = {arg[1].x,arg[1].y,arg[1].z,arg[1].w}
            end
        else
            v = {tonumber(arg[1])}
        end
    else
        v = arg
    end
    --local u = {}
    local n = 0
    for k,l in ipairs(v) do
        table.insert(self,l)
        --table.insert(u,l)
        n = n + 1
    end
    --self.vec = u
    self.size = n
    -- Shortcuts for low dimension
    self.x = self[1]
    self.y = self[2]
    self.z = self[3]
    self.w = self[4]
end

--[[
Test for zero vector.
--]]

function Vector:is_zero()
    for k,v in ipairs(self) do
        if v ~= 0 then
            return false
        end
    end
    return true
end

--[[
Test for equality.
--]]

function Vector:is_eq(u)
    if self.size ~= u.size then
        return false
    end
    for k,v in ipairs(self) do
        if v ~= u[k] then
            return false
        end
    end
    return true
end        

--[[
Inner product.
--]]

function Vector:dot(u)
    if self.size ~= u.size then
        return false
    end
    local d = 0
    for k,v in ipairs(self) do
        d = d +  v * u[k]
    end
    return d
end

--[[
Apply a given matrix (which is specified as a list of row vectors).
--]]

function Vector:applyMatrixLeft(m)
    if m.cols ~= self.size then
        return false
    end
    local u = {}
    for k,v in ipairs(m) do
        table.insert(u,self:dot(v))
    end
    return Vector(u)
end

function Vector:applyMatrixRight(m)
    if m.rows ~= self.size then
        return false
    end
    local u = {}
    local a = Vector.zero(m.cols)
    for k,v in ipairs(m) do
        a = a + self[k]*v
    end
    return a
end

--[[
Length of the vector
--]]

function Vector:len()
    local d = 0
    for k,v in ipairs(self) do
        d = d +  math.pow(v,2)
    end
    return math.sqrt(d)
end

--[[
Squared length of the vector.
--]]

function Vector:lenSqr()
    local d = 0
    for k,v in ipairs(self) do
        d = d +  math.pow(v,2)
    end
    return d
end

--[[
Norm infinity
--]]

function Vector:linfty()
    local m = 0
    for k,v in ipairs(self) do
        m = math.max(m,math.abs(v))
    end
    return m
end

--[[
Norm one
--]]

function Vector:lone()
    local m = 0
    for k,v in ipairs(self) do
        m = m + math.abs(v)
    end
    return m
end

--[[
Normalise the vector (if possible) to length 1.
--]]

function Vector:normalise()
    local l
    if self:is_zero() then
        print("Unable to normalise a zero-length vector")
        return false
    end
    l = 1/self:len()
    return self:scale(l)
end

--[[
Scale the vector.
--]]

function Vector:scale(l)
    local u = {}
    for k,v in ipairs(self) do
        table.insert(u,l*v)
    end
    return Vector(u)
end

--[[
Add vectors.
--]]

function Vector:add(u)
    if self.size ~= u.size then
        return false
    end
    local w = {}
    for k,v in ipairs(self) do
        table.insert(w, v + u[k])
    end
    return Vector(w)
end

--[[
Subtract vectors.
--]]

function Vector:subtract(u)
    if self.size ~= u.size then
        return false
    end
    local w = {}
    for k,v in ipairs(self) do
        table.insert(w, v - u[k])
    end
    return Vector(w)
end

--[[
Inline operators:

u + v
u - v
-u
u * v : scaling
u / v : scaling
u == v : equality
--]]

function Vector:__add(v)
    return self:add(v)
end

function Vector:__sub(v)
    return self:subtract(v)
end

function Vector:__unm()
    return self:scale(-1)
end

function Vector:__mul(v)
    if type(self) == "number" then
        return v:scale(self)
    elseif type(v) == "number" then
        return self:scale(v)
    else
        if self.is_a and self:is_a(Matrix) then
            return v:applyMatrixLeft(self)
        elseif v.is_a and v:is_a(Matrix) then
            return self:applyMatrixRight(v)
        end
    end
    return false
end

function Vector:__div(l)
    if type(l) == "number" then
        return self:scale(1/l)
    else
        return false
    end
end

function Vector:__eq(v)
    return self:is_eq(v)
end

function Vector:__concat(v)
    if type(v) == "table" 
        and v:is_a(Vector) then
            return self .. v:tostring()
        else
            return self:tostring() .. v
        end
end

function Vector:tostring()
    local t = {}
    for k,v in ipairs(self) do
        if type(v) == "table" and v.tostring then
            table.insert(t,v:tostring())
        else
            table.insert(t,v)
        end
    end
    return "(" .. table.concat(t,",") .. ")"
end

function Vector:__tostring()
    return self:tostring()
end

function Vector:tovec()
    if self.size == 2 then
        return vec2(self[1],self[2])
    elseif self.size == 3 then
        return vec3(self[1],self[2],self[3])
    else
        return vec4(self[1],self[2],self[3],self[4])
    end
end
        

function Vector.zero(n)
    if type(n) == "table" then
        n = n.size
    end
    u = {}
    for i = 1,n do
        table.insert(u,0)
    end
    return Vector(u)
end

function Vector:bitwiseReorder(b)
    local v
    if b then
        v = self
    else
        v = Vector({})
    end
    local h = BinLength(self.size)
    local m = math.pow(2,h)
    local l
    local w = {}
    for k = 1,m do
        if not w[k] then
            l = BinReverse(k-1,h)+1
            v[k],v[l] = self[l] or 0,self[k] or 0
            w[k],w[l] = 1,1
        end
    end
    v.size = m
    return v
end

function Vector:FFT(t)
    t = t or {}
    local v
    v = self:bitwiseReorder(t.inplace)
    local pi = math.pi
    if not t.inverse then
        pi = - pi
    end
    local r = v.size
    local n = BinLength(r) -1
    local i,j,d,s,m,fi,p,pr
    for k = 0,n do
        i = math.pow(2,k)
        j = math.pow(2,k+1)
        d = pi/i
        s = math.sin(d/2)
        m = Complex(-2 * s * s, math.sin(d))
        fi = Complex(1,0)
        for l = 1,i do
            for ll = l,r-1,j do
                p = ll + i
                pr = fi * v[p]
                v[p] = v[ll] - pr
                v[ll] = v[ll] + pr
            end
            fi = m*fi + fi
        end
    end
    return v
end
        

if cmodule.loaded "TestSuite" then
    
    testsuite.addTest({
        name = "Vector",
        setup = function()
            local w = Vector({1,2,3,4,5,6,7,8})
            local z = Vector({Complex(1,1), Complex(2,3), Complex(4,5)})
            for k,v in ipairs({
                {"Vector", w},
                {"Complex Vector", z},
                {"Bitwise Reverse", w:bitwiseReorder()},
                {"Bitwise Reverse", z:bitwiseReorder()},
                {"FFT", w:FFT()},
                {"FFT", z:FFT()}
            }) do
                print(v[1] .. ": " .. v[2])
            end
        end,
        draw = function()
        end
    })
end

return Vector

--]==]