Symbolic Math in Lua
Here's the Computer Algebra System in Lua I've started on.

Reference

Numbers

For the most part Lua numbers will work, and will automatically be replaced by symmath Constants (found in symmath.Constant). This is because Constant is a subclass of Expression, and I was too lazy to implement all the ridiculous number of edge cases required to handle Lua numbers, though I am thinking about doing this. Constant can be constructed with numbers (Constant(0), Constant(1), etc) or with complex values (Constant{re=1, im=0}, etc).

Complex values can be represented by symmath.complex, which uses builtin complex types when run within LuaJIT, or uses a pure-Lua alternative otherwise.

Some day I will properly represent sets (naturals, integers, rationals, reals, complex, etc), and maybe groups, rings, etc, but not for now.

Some day I will also add support for infinite precision or big integers, but not right now.

Variables

var = symmath.Variable(name[, dependencies])
var = symmath.var(name[, dependencies])
var1, var2, ... = symmath.vars(name1, name2, ...)
Create a variable with given name, and optionally a list of which variables it is dependent on for differentiation. By default variables of different names have a derivative of zero.

var:depends(var1, var2, ...)
Specify the variables that var is dependent on for differentiation.

func, code = symmath.compile(expr, {var1, var2, ...}, language)
func, code = expr:compile{var1, var2, ...}
Compiles an expression to a Lua function with the listed vars as parameters.
language can be one of the following:

Arithmetic

symmath.op.unm(a)
-a when used with symmath expressions
Creates an expression representing the negative of the parameter.

symmath.op.add(a, b)
a + b when used with symmath expressions
Creates an expression representing the sum of the parameters.

symmath.op.mul(a, b)
a * b when used with symmath expressions
Creates an expression representing the product of the parameters.

symmath.op.sub(a, b)
a - b when used with symmath expressions
Creates an expression representing the difference of the parameters.

symmath.op.div(a, b)
a / b when used with symmath expressions
Creates an expression representing a fraction of the parameters.

symmath.op.pow(a, b)
a ^ b when used with symmath expressions
Creates an expression representing the first parameter raised to the power of the second parameter.

symmath.op.mod(a, b)
a % b when used with symmath expressions
Creates an expression representing the first parameter modulo the second.

Equations

eqn = symmath.op.eq(lhs, rhs)
eqn = lhs:eq(rhs)
Creates an equation representing the equality lhs = rhs.

eqn = symmath.op.ne(lhs, rhs)
eqn = lhs:ne(rhs)
Creates an equation representing the inequality lhs ≠ rhs.

eqn = symmath.op.lt(lhs, rhs)
eqn = lhs:lt(rhs)
Creates an equation representing the inequality lhs < rhs.

eqn = symmath.op.le(lhs, rhs)
eqn = lhs:le(rhs)
Creates an equation representing the inequality lhs ≤ rhs.

eqn = symmath.op.gt(lhs, rhs)
eqn = lhs:gt(rhs)
Creates an equation representing the inequality lhs > rhs.

eqn = symmath.op.ge(lhs, rhs)
eqn = lhs:ge(rhs)
Creates an equation representing the inequality lhs ≥ rhs.

lhs = eqn:lhs()
Returns the left hand side of the equation.

rhs = eqn:rhs()
Returns the right hand side of the equation.

soln1, soln2, ... = symmath.solve(eqn, var)
soln1, soln2, ... = eqn:solve(var)
Returns solutions to the equation. If eqn is not an Equation then returns solutions to eqn = 0.

newexpr = expr:subst(eqn)
Shorthand for expr:replace(eqn:lhs(), eqn:rhs()).

Expressions

symmath.simplify(expr)
expr:simplify()
expr()
Simplifies the expression.

symmath.clone(expr)
Clones an expression.
This also replaces Lua numbers with symmath Constant objects.

symmath.replace(expr, find, repl, callback)
expr:replace(find, repl, callback)
Replaces portions of an expression with another.
expr = expression to change.
find = sub-expression to find.
repl = sub-expression to replace.
callback(node) = callback per node, returns 'true' if we don't want to find/replace this tree.

symmath.map(expr, callback)
expr:map(callback)
Maps sub-expressions in an expression to new sub-expressions.
expr = expression.
callback(node) = callback that returns nil if it leaves the tree untouched, returns a value if it wishes to change the tree.

symmath.eval(expr, {[var1]=value, var2name=value, ...})
expr:eval{[var1]=value, var2name=value, ...}
Calculates the numeric value of the expression.

symmath.polyCoeffs(expr, var)
expr:polyCoeffs(var)
Returns a table of coefficients with keys 0 through the degree of the polynomial, and 'extra' containing all non-polynomial terms.

Calculus

symmath.diff(expr, var1, var2, ...)
expr:diff(var1, var2, ...)
Differentiates the expression with respect to the given variable.

Linear Algebra

A = symmath.Matrix({expr11, expr12 ...}, {expr21, expr22, ...}, ...)
Create a matrix of expressions.

A = symmath.Array(...)
Create an array of expressions. Same deal as Matrix but with any arbitrary nesting depth, and without Matrix-specific operations.

AInv, I, message = A:inverse([b, callback, allowRectangular])
returns
AInv = the inverse of matrix A.
I = what's left of the Gauss-Jordan solver (typically the identity matrix).
message = any error that occured during solving.
args:
b = solution vector. If specified, returns x such that A x = b.
callback = function to be called after each Gauss-Jordan operation.
allowRectangular = set to true to allow inverting rectangular matrices.

APlus, determinable = A:pseudoInverse()
Returns the pseudoinverse of A (stored in APlus) and whether the pseudoinverse is determinable.

d = A:determinant()
Returns the determinant of A.

At = A:transpose()
Returns the transpose of A

I = Matrix.identity([m, n])
Returns 1 if no arguments are provided.
Returns a m x m identity matrix if one argument is provided.
Returns a m x n identity matrix if two arguments are provided.

D = Matrix.diagonal(d1, d2, d3, ...)
Returns a n x n matrix with diagonal elements set to d1 ... dn, for n the number of arguments provided.

tr = A:trace()
Returns the trace of A.

Tensors

Tensor.coords{ {variables={t,x,y,z}} }
Specifies that tensors will be using coordinates t,x,y,z

Tensor.coords{ {variables={t,x,y,z}, meric=g} }
Specifies that tensors will be using coordinates t,x,y,z with metric 'g' (a Matrix or 2D array). The metric inverse will be automatically computed.

Tensor.coords{ {variables={t,x,y,z}, meric=g, metricInverse=gU} }
Specifies that tensors will be using coordinates t,x,y,z with metric 'g' (a Matrix or 2D array) and metric inverse 'gU'.

Tensor.coords{ {variables={t,x,y,z}}, {symbols='ijklmn', variables={x,y,z}} }
Specifies that tensors will be using coordinates t,x,y,z, except for indexes ijklmn which will only use x,y,z.
At the moment conversion between maps is very ugly/incomplete.

Tensor.metric(g, [gU, symbol])
Specifies to use metric 'g' for the default coordinate system (assuming one has been defined with Tensor.coords).

t = Tensor'_abc'
Creates a degree-3 covariant tensor 't' with all values initialized to 0.

... tensor summation / multiplication ...
( u'^a' * v'_b' * eps'_a^bc_d' )() produces a degree-2 tensor perpendicular to 'u' and 'v' (assuming eps is defined as the Levi-Civita tensor).

... comma derivatives (semicolon derivatives almost there, just need to store the connection coefficients) ...

... assignment: S = T'_ab' - T'_ba'
in this case the indexes of 'S' are picked on a first-come, first-serve basis. If you want to be certain of the assignment, use the following:

S = Tensor'_ab'
S['_ab'] = (T'_ab' - T'_ba')()
(the final () is shorthand for :simplify(), which will evaluate the expression into the tensor structure).

... index gymnastics (so long as you defined a metric): v = Tensor('_a', ...) print(v'^a'()) will show you the contents of v^a = g^ab v_b.

t:print't'
Prints the tensor's contents.

t:printElem't'
Prints the individual nonzero values of the tensor, or '0' if they are all zero.

TODO

Download the Source Code Here
GitHub: http://github.com/thenumbernine/symmath-lua