Introduction  ::   complex  ::   matrix  ::   fft  ::   stat  ::   rng  ::   mathx

Introduction

What Numeric Lua is

Numeric Lua is a numerical package for the Lua programming language. It includes support for complex numbers, multidimensional matrices, random number generation, fast Fourier transforms, and special functions. Most of the routines are simple wrappers for well known numerical libraries: complex numbers and part of the extended math modules come from C99; other special functions, including statistical functions, are adapted from Netlib's SLATEC and DCDFLIB; random number generation is based on Takuji Nishimura and Makoto Matsumoto's Mersenne Twister generator as the "engine" (uniform deviates) and Netlib's RANLIB for the remaining deviates; fast Fourier transforms are implemented from FFTW; and the matrix package draws most of its numeric intensive routines from Netlib's ubiquitous BLAS and LAPACK packages.

Numeric Lua tries to maintain Lua's minimalist approach by providing bare-bone wrappers to the numerical routines. The user can use the outputs for further computations and is then fully responsible for the results. Other Lua features are also available, such as OO simulation through metamethods and functional facilities. A basic API is provided in order to promote extensibility. Also, check numlua.seeall for a quick way to start using Numeric Lua.

Numeric Lua is licensed under the same license as Lua -- the MIT license -- and so can be freely used for academic and commercial purposes.

numlua.buffer(opt [, arg])

Operates on Numlua buffer according to opt:

numlua.type(obj)

If obj is an userdata with metatable mt and mt.__type is not nil, returns mt.__type, otherwise returns the type of obj (type(obj)).

numlua.opmode ([mode])

Returns and optionally sets the operation mode to in-place. If mode is absent, returns the current operation mode; otherwise, sets the operation mode to mode and returns the old mode.

Note that metamethods are affected by opmode; for example, if numlua.opmode() is true,-a updates a in-place and the statement u = u + c * v updates both v and u.

numlua.seeall

A require call to numlua.seeall loads most of the functions in Numlua to the global table to make it easier to prototype and test code or even just use Numlua is a flexible computational shell. In particular, numlua.seeall does the following:

complex

Complex numbers

The complex module implements complex numbers according to the C99 standard. To create a complex number a + b * i, simply issue complex(a, b) or a + b * complex.i, or even simply a + b * i when using numlua.seeall. You can set the real and imaginary parts of complex z with z.r = x and z.i = x respectively. Since complex numbers are userdata, regular assignment is a reference assignment (z = c makes z a reference to c); to actually assign values, use z._ = c.

Similarly to the math library, all mathematical functions for C99 complex numbers are available; please check your man pages for more details about the following functions:

     abs      arg     acos    acosh     asin   asinh
    atan    atanh     conj      cos     cosh     exp
    imag      log   logabs     proj     real     sin
    sinh     sqrt      tan     tanh

In particular, to recover the real and imaginary parts of complex z, use z:real() and z:imag() respectively. In addition to these methods, add, sub, mul, div, and pow implement arithmetical methods and metamethods, that is, a + b is equivalent to complex.add(a, b) and so on. The "length" of z, #z, is complex.abs(z).

All methods accept an optional last argument -- whose default value is given by numlua.opmode -- that codes for in-place applications. For instance,

    > c = complex(2, -1)
    > print(c:exp(), c)
    3.9923240484413-6.217676312368i 2-1i
    > print(c:exp(true), c) -- in-place
    3.9923240484413-6.217676312368i 3.9923240484413-6.217676312368i

and while c = complex.exp(z * i) * a + b; print(c) creates (and discards) four temporary complex numbers, 

    op = numlua.opmode(true) -- in-place by default now
    c = complex(0, 1); c = complex.exp(c * z) * a + b
    numlua.opmode(op) -- restore previous mode

creates no temporary userdata as all operations are in-place. Note that operand order is important even when the operation is commutative: if numlua.opmode() is true, a + b is equivalent to complex.add(a, b, true) where the result is stored at a, while b + a, or complex.add(b, a, true), stores at b.

matrix

Multidimensional matrices

The matrix module implements multidimensional matrices. New matrix objects are created by specifying dimensions with matrix.new (uninitialized), matrix.zeros or matrix.ones (filled with zeros or ones), or by directly specifying entries in a table and using matrix (or matrix.fromtable.) Matrices can be real or complex, which can be set in the last argument of these functions. To retrieve the dimensions of a matrix use matrix.shape or matrix.size; in particular, m:size"#" returns the number of dimensions of a matrix m.

To access and set entries in a matrix m with n dimensions you can simply use m[i1]...[in]; for example, to set the entry at row i and column j of a two-dimensional matrix you can use m[i][j] = x. Alternatively, you can use matrix.get and matrix.set. These two functions can also be used to produce a reference to a submatrix and to assign whole submatrices; as an example, m:set(a) assigns a to m, either element-wise if a is a conformable matrix or by filling m with a if a is a number or complex. You can also use matrix.copy to copy a matrix. More general matrix references can be created with matrix.section or matrix.slice.

Entries in a matrix can be set more systematically using functional facilities: matrix.apply sets entries based on a function that takes indices as arguments while matrix.map takes a function on entries. Both functions operate in-place and in element-order, that is, in column-major order (first indices vary faster.) Element-order is also respected when iterating on a matrix with matrix.entries.

Logical methods -- matrix.find, matrix.ifelse, and matrix.which -- are helpful for filtering matrix entries. These functions always take a condition argument, which can be: a boolean function on each entry; or, if a number x, the boolean function function (e) return e == x end; or, if a conformable matrix x, function (e) return e == ex end where ex is the entry in x at the same element-order index of e, and so, a boolean function for element-wise comparison. The condition argument defaults to function (e) return e ~= 0 end.

Common mathematical operations can be carried using matrix.map but for convenience matrix provides the following functions:

      abs    acos    acosh    asin    asinh    atan
    atanh     cos     cosh     exp      log     sin
     sinh    sqrt      tan    tanh

All these functions take an optional second argument -- whose default value is given by numlua.opmode -- that specifies if the operation should be in-place. Thus, calling m:f() where f is any of the functions above, is semantically equivalent to matrix.map(inplace and m or m:copy(), f).

Arithmetical methods matrix.add, matrix.mul, matrix.div, and matrix.pow operate element-wise and accept an optional last argument that signals in-place application to the first argument. Matrix multiplication is implemented in three specialized methods: matrix.mmul for general matrix multiplication, matrix.hemul for Hermitian matrix multiplication, and matrix.trmul for triangular matrix multiplication. Arithmetical metamethods are, not surprisingly, based on these methods and should work as expected: __add and __sub are based on matrix.add, __mul is based on matrix.mul (for scalars) and matrix.mmul, __mod is matrix.ls (left division), and __div is based on matrix.div or matrix.ls for right division (a / b is equivalent to t(t(b) % t(a)) when a and b are conformable matrices, where t is matrix.t, the transpose operator.)

Note however, that even though metamethods make the code cleaner and more intuitive, it is often computationally more efficient to use the basic methods to avoid creation of temporary matrices. Here is a quick example:

    -- updates invA to inv(A + u * v') according to the
    -- Sherman-Morrison formula: inv(A + u * v') = inv(A)
    --          - (inv(A) * u * v' * inv(A)) / (1 + v' * inv(A) * u)
    function shermanmorrison (invA, u, v)
      local x = v * invA -- v' * inv(A)
      return invA:copy():mmul(invA * u, x, "n", "n", -1 / (1 + dot(x, u)))
    end

The next function is equivalent,

    function naiveshermanmorrison (invA, u, v)
      local x = v * invA -- v' * inv(A)
      return invA - invA * u * x / (1 + dot(x, u))
    end

but creates two extra temporary matrices: one from t = (invA * u) * x and other from t / (1 + dot(x,u)).

While favoring in-place operations can help avoid temporary object creation, the adoption of numlua.opmode can make notation lighter and less functional. It should be used with care, however, since operations being applied to the first argument of the methods break operation commutativity; for instance, while both matrix.add(a, b) and matrix.add(b, a) return a new matrix c = a + b, matrix.add(a, b, true) stores -- or better, applies -- the result in a and matrix.add(b, a, true) stores the result in b. Another example:

    function sweep (a, k)
      local ck = a:col(k)
      local b = ck:copy()
      local ak, d = a[k], b[k]
      ak:div(d, false, true) -- a[k] = a[k] / b[k], in-place
      for i = 1, #a do
        if i ~= k then
          a[i]:add(ak, -b[i], true) -- a[i] = a[i] - b[i] * a[k], in-place
        end
      end
      ck._ = b:div(-d, false, true) -- ck = b, b = -b / d in-place
      ck[k] = 1 / d
      return a
    end

and now using in-place operations by default and metamethods:

    numlua.opmode(true) -- in-place
    function sweep (a, k)
      local ck = a:col(k)
      local b = ck:copy()
      local ak, d = a[k], b[k]
      ak = ak / d
      for i = 1, #a do
        if i ~= k then
          a[i]:add(ak, -b[i]) -- a[i] = a[i] - ak * b[i]
        end
      end
      ck._ = b / (-d) -- equivalent to ck:set(b / (-d))
      ck[k] = 1 / d
      return a
    end

Finally, matrices can be serialized (in HDF5 format) using matrix.save and retrieved using matrix.load.

matrix.get(m, i1 [, i2 [... [, in]]])
matrix.get(m, eoindex)

Returns a submatrix of m. There are two ways of specifying the submatrix:

See also: matrix.set

matrix.set(m, i1 [, i2 [... [, in]]], value)
matrix.set(m, eoindex, value)
matrix.set(m, what, value)

Sets a submatrix of m to value. There are three ways of specifying the submatrix:

Parameter value can be a number (or complex) or a conformable matrix according to the submatrix being indexed.

See also: matrix.get

matrix.slice(m [, first] [, last] [, step])

Returns a section of matrix m that is a slice in the first dimension of m starting at first (defaults to 1), ending at last (defaults to #m), and alternating by step (defaults to 1) positions. The returned section is a reference to m.

A shortcut for matrix.slice(m, first, last, step) is m(first, last, step).

See also: matrix.section

matrix.section(m, section)

Returns a section of matrix m according to table section. section is taken as a table with integer indices from 1 to n, the number of dimensions of m; the i-th value of section should contain another table specifying a slice in the form {first, last, step}, where first and step default to 1 and last defaults to m:size(i). The returned section is a reference to m.

A shortcut for matrix.section(m, section) is m(section).

Example:

    -- assume a and b are real two-dimensional matrices
    -- check `matrix.kronecker` for a more efficient version
    function kronecker (a, b)
      local ra, ca = a:shape()
      local rb, cb = b:shape()
      local x = matrix.new(ra * rb, ca * cb)
      for i = 1, ra do
        for j = 1, ca do
          x{{(i-1)*rb+1, i*rb}, {(j-1)*cb+1, j*cb}}._ = b * a[i][j]
        end
      end
      return x
    end

See also: matrix.slice

matrix.add(a, b [, s [, inplace]])

Adds s * b to a, where s defaults to 1. If b is a number s is ignored and sets a[i] = a[i] + real(b) if a is real, or a[i] = a[i] + b if a is complex; otherwise, sets a[i] = a[i] + b[i] * s if b is a conformable matrix. In particular, matrix.add(a, b, -1) subtracts b from a.

If inplace is true, a is updated in-place.

matrix.mul(a, b [, inplace])

Performs element-wise multiplication of a and b. If b is a number, sets a[i] = a[i] * real(b) if a is real, or a[i] = a[i] * b if a is complex; otherwise, sets a[i] = a[i] * b[i] if b is a conformable matrix.

If inplace is true, a is updated in-place.

matrix.div(a, b [, right [, inplace]])

Performs element-wise division of a and b. If b is a number, sets a[i] = a[i] / real(b) if a is real, or a[i] = a[i] / b if a is complex; otherwise, sets a[i] = a[i] / b[i] if b is a conformable matrix.

If right is true, a is divided in the right, that is, a = b / a. If inplace is true, a is updated in-place.

matrix.pow(a, b [, inplace])

Performs element-wise power of a and b. If b is a number, sets a[i] = a[i] ^ real(b) if a is real, or a[i] = a[i] ^ b if a is complex; otherwise, sets a[i] = a[i] ^ b[i] if b is a conformable matrix.

If inplace is true, a is updated in-place.

matrix.new(d1, d2, ... [, iscomplex])
matrix.zeros(d1, d2, ... [, iscomplex])
matrix.ones(d1, d2, ... [, iscomplex])

Returns a matrix with dimensions d1, d2, and so on. If last argument iscomplex is true, the returned matrix is complex. The entries are arbitrary (not initialized) when using new, initialized to zero when using zeros, and initialized to one when using ones.

matrix.eorder(m, i1, i2, ..., in)

Returns the element-order index in m that corresponds to the dimensional indexes i1, i2, ..., in, where n is the dimensionality of m.

See also: matrix.eindex

matrix.eindex(m, eopos)

Returns the dimensional indexes in m that correspond to element-order index eopos.

Example:

    function list (m)
      for i, e in m:entries(true) do
        local t = {m:eindex(i)}
        t[#t + 1] = e
        print(unpack(t))
      end
    end

See also: matrix.eorder, matrix.entries

matrix.entries(m [, eorder])

Returns an iterator over the entries of m so that the following construction iterates over the entries of m:

    for i1, i2, ..., in, e in m:entries() do
      -- body
    end

Here n = m:size"#" is the number of dimensions of m and e = m[i1][i2]...[in] takes the value of each entry in m.

This construction is semantically equivalent to

    for i = 1, m:size"*" do
      local i1, i2, ..., in = m:eindex(i)
      local e = m:get(i1, i2, ..., in)
      -- body
    end

See also: matrix.eindex

matrix.size(m [, index])

Returns the size of m, depending on parameter index:

The default value for index is 1.

See also: matrix.shape

matrix.shape(m [, index [, complex]])

Returns the dimensions of m starting at dimension index (defaults to 1) and if m is complex when complex is true.

See also: matrix.size

matrix.iscomplex(m)

Returns true if m is complex and false otherwise.

matrix.real(a)

Returns a real matrix that references the real part of a if a is complex or returns a itself if a is real.

See also: matrix.imag, matrix.complex

matrix.imag(a)

Returns a real matrix that references the imaginary part of a if a is complex or returns matrix.zeros(a:shape()) if a is real.

See also: matrix.real, matrix.complex

matrix.complex(a [, b])

Returns a complex matrix m such that real(m) = real(a) and imag(m) = imag(b). The default value for b is a, so, in particular, complex(a) returns a copy of a if a is complex. An error is raised if a and b are not conformable.

See also: matrix.real, matrix.imag

matrix.conj(a [, inplace])

Sets each element in a to its conjugate in-place if inplace is true, or returns the conjugate of a. If a is not complex, does nothing.

matrix.copy(m [, x])

Returns a copy of m and optionally sets its content to x, a number. matrix.copy is semantically equivalent to

    function matrix.copy (m, x)
      local c = matrix.new(matrix.shape(m, 1, true)) -- frame
      return set(c, x or m)
    end

See also: matrix.set

matrix.reshape(m, d1, d2, ..., dn)

Returns a reference to m that has dimensions d1, d2, ..., dn. The new dimensions need to be consistent, that is, d1 * d2 * ... * dn == m:size"*".

See also: matrix.size

matrix.spread(m [, dim [, count]])

Returns a replicated version of m along dimension d (defaults to 1) count number of times (defaults to 1).

matrix.find(m [, cond [, reverse]])

Returns the first element-order index in m whose respective entry satisfies cond. If reverse is true, returns the first occurrence in m that does not satisfy cond. If no entry in m satisfies cond, nil is returned. Note that matrix m is traversed only up to the first successful entry.

Examples:

    function any (m, cond)
      return matrix.find(m, cond) ~= nil
    end

    function all (m, cond)
      return matrix.find(m, cond, true) == nil
    end

matrix.ifelse(m, cond, x [, y])

Sets the entries in m that satisfy cond to x and, if y is specified, the entries in m that do not satisfy cond to y. Parameters x and y can be numbers or conformable matrices to m.

Examples:

    -- note that matrix.which(m, cond, "#") is preferred
    function count (m, cond)
      return matrix.sum(matrix.ifelse(m:copy(), cond, 1, 0))
    end

    -- mask has ones for true and zeros for false
    function merge (x, y, mask)
      return matrix.ifelse(mask:copy(), 1, x, y)
    end

    local lt = function (x, y) return x < y and 1 or 0 end
    function min (x, y)
      return ifelse(matrix.map(x:copy(), y, lt), 1, x, y)
    end
    function max (x, y)
      return ifelse(matrix.map(x:copy(), y, lt), 0, x, y)
    end

See also: matrix.which

matrix.which(m [, cond [, what]])

Returns entries in m that satisfy cond depending on what:

Parameter what defaults to "k".

Examples:

    -- simple alias
    function count (m, cond)
      return matrix.which(m, cond, "#")
    end

    -- mask has ones for true and zeros for false
    function pack (m, mask)
      return matrix.which(m, mask, "v")
    end
    function unpack (v, mask, m)
      return matrix.set(m, matrix.which(m, mask), v)
    end

See also: matrix.ifelse

matrix.apply(m, f [, eorder])

Applies a function f to each entry in m in-place, that is, sets m[i1,...,in] = f(i1, ..., in, m[i1, ..., in]) where i1, ..., in are indices in m (m has dimension n). If eorder is true, f takes indices in element-order. If f does not return a numeric value, m is not changed for the called index.

Examples:

    local h = function (i, j) return 1 / (i + j - 1) end
    function hilbert (n)
      return matrix.new(n, n):apply(h)
    end

    local circ = function (v, n)
      return function (i, j) return v[(j - i) % n + 1] end
    end
    function circulant (v)
      local n = #v
      return matrix.new(n, n):apply(circ(v, n))
    end

matrix.fold(m, f [, init])

Computes the left fold of function f over the entries in m in element-order, that is, returns

    f(...f(f(init, m[1]), m[2])..., m[n])

See also: matrix.sum

matrix.map(m1, m2, ..., mn, f)

Sets, in element-order, the i-th entry in m1 in-place to

    m1[i] = f(m1[i], m2[i], ..., mn[i])

if the returned value is a number or complex, where m2, ..., mn are all conformable to m1.

matrix.sum(m [, alpha [, init]])

Performs a linear fold on matrix m in element-order, that is, computes

    (...((init * alpha + m[1]) * alpha + m[2]) ... ) * alpha + m[n])

where n = m:size"*".

Parameters alpha and init default to 1 and 0 respectively; in particular, matrix.sum(m) computes the simple sum of the entries in m. It is semantically equivalent to

    matrix.fold(m, function(a, e) return a * alpha + e end, init)

See also: matrix.fold

matrix.min(m)

Returns the minimum element in m and the minimizing element-order index of m. Complex numbers are compared based on their real parts, and their imaginary parts are only used to break ties.

See also: matrix.max, matrix.sort

matrix.max(m)

Returns the maximum element in m and the maximizing element-order index of m. Complex numbers are compared based on their real parts, and their imaginary parts are only used to break ties.

See also: matrix.min, matrix.sort

matrix.sort (m [, decreasing [, returnindex]])

Sorts m in-place. If decreasing is true, sorts in decreasing order. If returnindex is true returns index such that m[i] was m[index[i]] before sorting, where i is in element-order; otherwise, returns m.

As in matrix.min and matrix.max, complex numbers are compared based on their real parts, and their imaginary parts are only used to break ties.

Example:

    iv = matrix.sort(v, false, true) -- increasing, return index
    x = matrix.new(matrix.shape(v, 1, true))
    x[iv] = v -- restore v

See also: matrix.min, matrix.max

matrix.linspace(a, b [, n])

Returns a vector v of length n -- taken to be int(|b - a| + 1) if not specified -- such that v[i] = (i - 1) * s where s is (b - a) / (n - 1), and so v[1] = a and v[n] = b.

Examples:

    v = linspace(1, n) -- v = 1:n

    function arith (a, r, n)
      return linspace(a, r * (n - 1) + a, n)
    end

    function seq (a, b, step)
      local n = floor(abs((b - a) / s + 1))
      return linspace(a, b, n)
    end

matrix.dot(a, b [, trans])

Computes the dot product of conformable matrices a and b, a^H b as the sum of their element-wise product. If trans is true, the conjugate of a is not used (and so a^T b is actually computed.)

matrix.cross(a, b [, inplace])

Computes the cross product of ternary vectors a and b in-place in a if inplace is true.

matrix.col(m, i)

Returns a reference to the i-th column of the two-dimensional matrix m.

matrix.transpose (m [, hermitian])

Returns the transpose of matrix m. If hermitian is true, returns the conjugate transpose of m.

matrix.diag(m [, shift])

If m is a vector, returns a square matrix such that the shift diagonal of the matrix is m; otherwise, if m is two-dimensional, returns a reference to the shift diagonal of m. If shift is zero (the default value), the main diagonal is referenced; if shift > 0, the shift diagonal above the main diagonal is referenced, and if shift < 0, the shift diagonal below the main diagonal is referenced.

Example:

    function trace (m)
      local r, c = m:shape()
      assert(m:size"#" == 2 and r == c, "square matrix expected")
      return m:diag():sum()
    end

matrix.concat(m1, ..., mn [, colcat])

Concatenates conformable matrices m1, ..., mn row-wise. If colcat is true, concatenates column-wise.

matrix.c(v1, ..., vn)

Concatenates vectors or numbers v1, ..., vn into a new vector.

matrix.swap (a, b)

Swaps the contents of conformable vectors a and b.

matrix.pivot(m, pvt [, colpivot [, inplace]])

Given pivot vector pvt, returns a new matrix p with the same rows of m but with rows i and pvt[i] swapped. Note that, in particular, p = P * m, where P = eye(#m):pivot(pvt), that is, P is the permutation matrix for pvt.

If colpivot is true, columns are swapped instead of rows, that is, returns p = m * P where P = eye(#m):pivot(pvt, true).

If inplace is true, rows/columns are swapped in-place.

Examples:

    function perm (p) -- column permutation
      local n = #p
      local s = matrix.linspace(1, n):pivot(p)
      local P = matrix.zeros(n, n)
      for i = 1, n do P[i][s[i]] = 1 end
      return P
    end

    -- see `matrix.lu` and `matrix.qr` for details: 
    require "numlua.seeall"
    l, u, p = lu(a)
    print(pretty(perm(p) * a)) -- = a:pivot(p)) = l * u
    q, r, p = qr(a,true)
    print(pretty(q * r * t(perm(p)))) -- = (q * r):pivot(p,true) = a

matrix.trmul (x, A [, uplo [, invert [, trans [, side]]]])

For conformable two-dimensional matrices A and x, A square, returns the result of tri(A) * x in-place in x, where tri(A) takes the lower triangle section of A if uplo == "l" or uplo == "L", or the upper triangle section of A if uplo == "u" or uplo == "U". The default value for uplo is "L".

If invert is true, multiplies by the inverse of A.

If trans == "t" or trans == "T", returns tri(A)^T * x; if trans == "c" or trans == "C", returns tri(A)^H * x (conjugate transpose); if trans == "n" or trans == "N", does nothing (default behavior).

If side == "l" or side == "L", A is multiplied at the left of x (default); if side == "r" or side == "R", A is multiplied at the right, that is, returns x * tri(A).

For examples, see matrix.lu and matrix.qr.

matrix.hemul (C, A, [, inner, [, what [, alpha]]])

Given a Hermitian matrix C, performs the Hermitian update C = C + alpha * A^H * A if inner is true, or C = C + alpha * A * A^H otherwise. The default value for alpha is 1.

Parameter what specifies how C is to be updated; if what == "l" or what == "L", only the lower triangle of C is updated; if what == "u" or what == "U", only the upper triangle of C is updated; if what == "f" or what == "F", the lower triangle of C is updated by the multiplication and reflected to the strict upper triangle. The default value for what is "F".

Example:

    -- generates Wishart deviates
    function wishart (S, n)
      local m = #S
      local L = assert(chol(S, "L")) -- Cholesky factor: S = L * L^T
      local c = zeros(m) -- cache
      return function ()
        local w = zeros(m, m)
        for i = 1, n - 1 do
          c = trmul(rnorm(0, 1, c), L) -- c ~ N(0, S)
          w = hemul(w, c, false, "L") -- w = w + c * c^T
        end
        return hemul(w, trmul(rnorm(0, 1, c), L)) -- n-th run: full hemul
      end
    end

matrix.mmul (C, A, B [, transA, transB, alpha])

For conformable matrices A, B, and C, updates C with the multiplication C = C + alpha * transA(A) * transB(B). Parameters transA and transB specify possible operations for A and B:

If B is a vector, transB is ignored. Moreover, A can only be a vector if B is also a vector, in which case the outer product A * B^H is computed regardless of transA and transB.

matrix.norm(m [, what])

Computes the norm of m according to parameter what:

The default value for what is "F". If what == "m", what == "M", what == "i", or what == "I", also returns the index the maximizes the norm.

matrix.chol(A [, what [, inplace]])

Computes the Cholesky factorization of symmetric positive-definite matrix A. Parameter what defines the computation:

If an internal error occurs, returns nil and an error message; if A is not positive-definite, returns false and a warning message. If inplace is true, overwrites A with its Cholesky factor.

Examples:

    function matrix.isposdef (m)
      local c, msg = chol(m)
      if c == nil then error(msg) end
      return not c == false
    end

    -- generates multivariable normal deviates (also check rng.rmvnorm)
    function mvnorm (mu, S)
      local L = assert(chol(S, "L")) -- Cholesky factor: S = L * L^T
      return function (dest)
        local dest = dest or zeros(#mu)
        local s = rnorm(0, 1, dest) -- s ~ N(0, I_n)
        s = trmul(s, L) -- s = L * s, s ~ N(0, S)
        s = add(s, mu, true) -- s = s + mu, s ~ N(mu, S)
        return s
      end
    end

See also: matrix.hemul, matrix.trmul

matrix.lu(A [, inplace])

Computes the LU decomposition of A: returns an unit lower triangular matrix L, an upper triangular matrix U and a permutation vector pvt such that A:pivot(pvt) = L * U, or equivalently, A = P^T * L * U where P = eye(#A):pivot(pvt). If inplace is true, both L and U are returned in-place at A, L on the strict lower triangle and U on the upper triangle.

Examples:

    -- compute determinant of m
    function det (m)
      local c = assert(lu(copy(m), true)) -- in-place on copy
      return prod(diag(c))
    end

    -- solve linear system A * x = b
    function linsolve (A, b)
      local l, u, p = lu(A)
      local x = b:pivot(p) -- x = eye(#A):pivot(p) * b
      x:trmul(l, 'l', true) -- x = inv(l) * x
      x:trmul(u, 'u', true) -- x = inv(u) * x
      return x
    end

See also: matrix.pivot

matrix.rcond(A [, what [, inplace]])

Estimates the reciprocal of the condition number of matrix A in 1-norm, depending on parameter what:

If inplace is true and what codes for A as being Hermitian positive-definite or general, A is overwritten with either a Cholesky or LU decomposition respectively.

See also: matrix.svd, matrix.inv

matrix.inv(A [, what [, inplace [, norcond]]])

Returns the inverse of square matrix A and an estimate of the reciprocal of its condition number, according to parameter what:

If inplace is true and what codes for A as being Hermitian positive-definite or general, A is overwritten with either a Cholesky or LU decomposition respectively.

If norcond is true, the reciprocal of the condition number is not estimated.

Example:

    -- Square-root of a matrix A by Denman-Beavers algorithm
    function dbsqrtm (A, tol, maxiters)
      local opmode = numlua.opmode(true) -- set in-place operations
      local iY, iZ = zeros(A:shape()), zeros(A:shape()) -- buffers
      local s, n = inf, 1 -- norm(Y), #iterations
      local Y, Z = copy(A), eye(#A) -- initialize
      while true do
        iY._, iZ._ = Y, Z
        Y, Z = (Y + inv(iZ)) / 2, (Z + inv(iY)) / 2 -- in-place
        -- check termination:
        local f = norm(Y)
        if abs(f - s) <= tol or n > maxiters then break end
        s, n = f, n + 1
      end
      numlua.opmode(opmode) -- restore previous mode
      return Y, Z, s
    end

See also: matrix.rcond

matrix.svd(A [, what])

Computes the singular value decomposition of m-by-n matrix A. The SVD is specified as A = U * S * V^H, where U and V are unitary matrices whose first min(m, n) columns represent the left and right singular vectors of A, respectively, and S is a diagonal matrix whose min(m, n) diagonal entries contain the singular values of A. U, V, and S are computed according to parameter what:

The singular values are returned in descending order. Note that V^H is returned, not V.

If an internal error occurs, returns nil and a diagnostic message; if the algorithm fails to converge, returns false and a warning message.

Examples:

    function cond (m)
      local s = assert(svd(m, "n")) -- just singular values
      return s[1] / s[#m] -- condition number
    end

    local lt = function (x) return function(e) return e < x end end
    function matrix.rank (m, tol) -- effective rank according to tol
      local s = assert(svd(m, "n")) -- just singular values
      local tol = tol or 0
      if tol <= 0 then -- set default tolerance?
        tol = max(m:shape()) * eps * s[1]
      end
      local r = s:find(lt(tol))
      if r == nil then -- all > tol?
        return #s -- full rank: min(m:shape())
      end
      return r - 1 -- last position > tol
    end

For other examples, see matrix.null, matrix.orth, and matrix.pinv in matrix.lua.

See also: matrix.rcond

matrix.qr(A [, permute [, inplace]])

Computes the QR factorization of matrix A, that is, returns an unitary matrix Q and an upper trapezoidal matrix R such that A = Q * R.

If permute is true, computes the QR factorization with column pivoting by also returning a permutation vector p such that A = (Q * R):pivot(p, true), that is, P = eye(A:size(2)):pivot(p, true) satisfies A * P^T = Q * R.

In inplace is true, A is overwritten with R.

Example:

    -- solve linear system A * x = b
    function linsolve (A, b)
      local q, r = qr(A)
      local x = zeros(#b):mmul(A, b, 't') -- x = t(A) * b
      x:trmul(r, 'u', true, 't') -- x = t(inv(r)) * x
      x:trmul(r, 'u', true) -- x = inv(r) * x
      return x
    end

matrix.eig(A [, what [, hermitian]])

Computes the eigenvalues of A, and optionally the eigenvectors according to parameter what:

Note that the right eigenvectors of A are the left eigenvalues of A^H and vice-versa. If hermitian is true, A is assumed to be Hermitian and only its upper triangle is accessed.

When computing the eigenvalues -- and optionally eigenvectors -- of A, the matrix is balanced to improve accuracy; see matrix.balance for a description of the operations. There is currently no way of preventing balancing.

Examples:

    -- condition number
    function cond (m)
      local x = zeros(#m, #m):hemul(m) -- x = m * t(m)
      local e = eig(x, "n", true) -- sqrt(e) are singular values of m
      return math.sqrt(math.abs(max(e) / min(e)))
    end

    -- if F[n] is the n-th Fibonacci number,
    -- {{F[n-1], F[n]}, {F[n], F[n+1]}} = {{0, 1}, {1, 1}} ^ n
    local s, V = eig(matrix{{0, 1}, {1, 1}}, "r", true)
    local v = V[1] * V[2] -- element-wise multiplication
    local w = zeros(2) -- workspace
    function fib (n)
      w._ = s; w:pow(n, true) -- w = s ^ n, in-place
      return dot(v, w) -- V * diag(s ^ n) * t(V)
    end

See also: matrix.balance, matrix.svd

matrix.balance (A [, what])

Balances a matrix A as a way to reduce its 1-norm and thus improve the accuracy of the computed eigenvalues and eigenvectors. The operations involve, first, permuting A by a similarity transformation to isolate eigenvalues in the first 1 to l and last h to #A elements of the diagonal, and then scaling by a diagonal similarity transformation to rows and columns l + 1 to h - 1 to make them as close in norm as possible. Parameter what controls which operations are performed to balance A:

See also: matrix.eig

matrix.ls (A, b [, svd [, tol [, inplace]]])

Computes the least-squares solution to the linear system A * x = b, that is, finds x that minimizes the 2-norm of b - A * x, where A is a m-by-n matrix. Parameter b can be a (m-by-nrhs) matrix containing many right-hand side vectors as columns.

If svd is true, the solution is found using a singular value decomposition; otherwise, a QR factorization is applied. The effective rank of A is relative to tol, which defaults to max(m, n) * mathx.eps (see Example in matrix.svd for details).

If inplace is true, the solution is returned in b if the system is not underdetermined, that is, if A has at least as many rows as columns. In this case, b:slice(1,n) contains the solutions x and the sum of squares for the i-th column of b:slice(n+1, m) is the residual sum of squares for the solution of the i-th right-hand side.

The solution and the effective rank are returned; in addition, if svd is true, the singular values of A are also returned (as a vector).

Example:

    -- basic LS linear model fitting
    function matrix.lm (a, b, svd)
      local m, n = a:shape()
      assert(m >= n, "system is underdetermined")
      assert(b:size"#" == 1, "single RHS expected")
      local x, rank = ls(a, b, svd)
      -- report summary statistics
      local coef = x:slice(1, n)
      local rss = (b - a * coef):norm() ^ 2
      local rss0 = (b - b:sum() / m):norm() ^ 2
      local df = m - rank
      local F = df / (rank - 1) * (rss0 / rss - 1)
      local pvalue = 1 - stat.pf(F, rank - 1, df)
      return {coef = coef, rss = rss, df = df, F = F, pvalue =  pvalue}
    end

See also: matrix.svd, matrix.qr

matrix.fft(m [, inverse [, wisdomonly [, inplace]]])

Computes the fast Fourier transform of complex matrix m. If inverse is true, computes the normalized inverse FFT. wisdomonly is a FFTW-specific parameter: if true, a plan to compute the transform is only created if wisdom is available for the given problem, that is, a FFTW_WISDOM_ONLY flag is passed when building the plan. The default flag is FFTW_ESTIMATE. If inplace is true, m is updated in-place.

See also: matrix.fct, fft.plan

matrix.fct(m [, inverse [, wisdomonly [, inplace]]])

Computes the discrete cosine transform of real matrix m. If inverse is true, computes the normalized inverse DCT. wisdomonly is a FFTW-specific parameter: if true, a plan to compute the transform is only created if wisdom is available for the given problem, that is, a FFTW_WISDOM_ONLY flag is passed when building the plan. The default flag is FFTW_ESTIMATE. If inplace is true, m is updated in-place.

See also: matrix.fft, fft.plan

matrix.save(m, filename)

Saves matrix m in HDF5 file filename.

matrix.load(filename)

Loads a matrix from HDF5 file filename and returns the loaded matrix.

fft

Fast Fourier transforms

The fft module implements fast Fourier transforms. The simple interface to FFTs comprehends just matrix.fft and matrix.fct which operate directly on a matrix. If, however, you intend to execute the same transform many times it might be more efficient to use a plan, so check fft.plan. If more control on how plans are executed is needed, there is also an interface to the wisdom of FFTW in fft.wisdom.

Example:

    > a = matrix({{1, 2, 3}, {0, -1, -4}}, true)
    > b = matrix(a:shape(1, true))
    > p = fft.plan(b, false, fft.flag["measure"]) -- direct transform
    > ip = fft.plan(b, true, fft.flag["measure"]) -- inverse transform
    > b._ = a
    > p(); print(b:pretty()) -- execute direct
                      1+0i    1-1.7320508075689i    1+1.7320508075689i
                     11+0i   -4+3.4641016151378i   -4-3.4641016151378i
    > ip(); print(b:pretty()) -- execute inverse
        1+0i    2+0i    3+0i
        0+0i   -1+0i   -4+0i

fft.plan(m [[, inverse], flag])

Returns a FFTW plan that uses matrix m as a buffer to apply discrete transforms in-place in m. If m is complex, a fast Fourier transform is executed, otherwise, if m is real, a discrete cosine transform is performed. If inverse is true, normalized inverse transforms are computed.

Parameter flag specifies a planner flag. Possible values are available in table fft.flag. The default flag is fft.flag["measure"], FFTW_MEASURE. Note that only "estimate" (FFTW_ESTIMATE) and "wisdomonly" (FFTW_WISDOM_ONLY) do not overwrite buffer m; for all other flags, m has to be initialized after the plan is created.

If plan = fft.plan(m, ...), #plan returns m and plan() executes the plan.

See also: matrix.fft, matrix.fct

fft.wisdom([wisdom])

Provides an interface for FFTW wisdoms. If wisdom is nil, returns the currently accumulated wisdom; if wisdom is true, forgets accumulated wisdom; if wisdom is a string, imports wisdom from it.

stat

Statistics and probability methods

The stat module implements statistical routines: probability mass and density functions, cumulative distribution functions (cdf) and inverse cdf (quantile) functions for a number of distributions; and factors.

stat.factor(t)

Constructs a factor based on table t, a sequence. A factor is a more compact representation of t where each distinct entry -- a level -- is assigned an integer key, the key level. If f is a factor built from t, f[i] returns t[i], f(i) returns the key level of t[i], and f() returns the levels of f as a table; in particular, f[i] is equivalent to f()[f(i)]. The expression #f returns a vector with level counts, that is, (#f)[l] returns the number of entries in t with key level l and thus values f()[l]; note that (#f):sum() is #t.

Example:

    > f = stat.factor{"a", "c", "a", "a", "b", "c"}
    > for l, level in ipairs(f()) do print(l, level, (#f)[l]) end
    1       a       3
    2       c       2
    3       b       1
    > for i = 1, (#f):sum() do print(i, f[i], f(i)) end
    1       a       1
    2       c       2
    3       a       1
    4       a       1
    5       b       3
    6       c       2

See also: factor.fold, factor.partition, factor.design

factor.fold(f, m, func [, init])

Folds the entries in matrix m according to the levels in factor f and fold function func, that is, for each key level l, returns a vector v such that v[l] = func(m[i1], func(m[i2], func(..., func(m[in], init)...))) where i1, ..., in satisfy f(i1) = ... = f(in) = l. The default value for init is zero.

Example:

    > f = stat.factor{"a", "c", "a", "a", "b", "c"}
    > add2 = function(x, y) return x + y end
    > x = matrix.linspace(1, (#f):sum())
    > for l, v in f:fold(x, add2):entries() do print(l, v) end
    1       8
    2       8
    3       5

See also: stat.factor, factor.partition

factor.partition(f, m)

Partitions the matrix m according to the levels in f by returning a table t such that t[l] is a vector containing the entries i1, ..., in in m such that f(i1) = ... = f(in) = l.

Example:

    > f = stat.factor{"a", "c", "a", "a", "b", "c"}
    > x = matrix.linspace(1, (#f):sum())
    > vprint = function(v) return table.concat(v:totable(), ", ") end
    > for l, v in ipairs(f:partition(x)) do print(l, vprint(v)) end
    1       1, 3, 4
    2       2, 6
    3       5

See also: stat.factor, factor.fold

factor.design(f, [ref])

Retuns design matrix for factor f. If ref is zero (default value) then no contrasts are assumed; otherwise, return a "sum to zero" contrast matrix with ref as reference level.

Example:

    > f = stat.factor{"a", "c", "a", "a", "b", "c"}
    > print(f:design():pretty())
       1   0   0
       0   1   0
       1   0   0
       1   0   0
       0   0   1
       0   1   0
    > print(f:design(2):pretty())
        1    1    0
       -1    1   -1
        1    1    0
        1    1    0
        0    1    1
       -1    1   -1

stat.dbeta(x, a, b)
stat.pbeta(x, a, b)
stat.qbeta(p, a, b)

Computes the density (dbeta), cumulative distribution function (pbeta), and inverse cdf (qbeta) for the beta distribution with parameters a >= 0 and b >= 0. The density is given by x^(a - 1) * (1 - x)^(b - 1) / B(a, b) for 0 <= x <= 1, where B(a, b) is the beta function.

See also: mathx.beta

stat.dbinom(s, xn, pr)
stat.pbinom(s, xn, pr)
stat.qbinom(p, xn, pr)

Computes the probability mass function (dbinom), cumulative distribution function (pbinom), and inverse cdf (qbinom) for the binomial distribution with size xn >= 0 and probability of success on each trial 0 <= pr <= 1. The density is given by choose(xn, s) * pr^s * (1 - pr)^(xn - s) for s = 0, ..., xn, where choose(xn, s) is the binomial coefficient.

stat.dchisq(x, df [, pnonc])
stat.pchisq(x, df [, pnonc])
stat.qchisq(p, df [, pnonc])

Computes the density (dchisq), cumulative distribution function (pchisq), and inverse cdf (qchisq) for the chi-squared distribution with df > 0 degrees of freedom and non-centrality parameter pnonc. If pnonc ~= 0, df can be zero.

stat.dexp(x [, r])
stat.pexp(x [, r])
stat.qexp(p [, r])

Computes the density (dexp), cumulative distribution function (pexp), and inverse cdf (qexp) for the exponential distribution with rate r. The default value for r is one. The exponential distribution has density r * exp(-r * x).

stat.df(x, dfn, dfd)
stat.pf(x, dfn, dfd [, pnonc])
stat.qf(p, dfn, dfd [, pnonc])

Computes the density (df), cumulative distribution function (pf), and inverse cdf (qf) for the F distribution with dfn >= 0 and dfd >= 0 degrees of freedom in the numerator and denominator, and non-centrality parameter pnonc.

stat.dgamma(x, shape [, rate])
stat.pgamma(x, shape [, rate])
stat.qgamma(p, shape [, rate])

Computes the density (dgamma), cumulative distribution function (pgamma), and inverse cdf (qgamma) for the gamma distribution with parameters shape >= 0 and rate >= 0. The default value for rate is one. A gamma distribution with shape s and rate r has density r^s / G(s) * x^(s - 1) * exp(-r * x) where G is the gamma function, for x >= 0.

See also: mathx.gamma

stat.dhyper(x, r, b, n)
stat.phyper(x, r, b, n)

Computes the probability mass function (dhyper) and cumulative distribution function (phyper) for the hypergeometric distribution with parameters r, b >= 0 and 0 <= n <= r + b. The density is given by choose(r, x) * choose(b, n - x) / choose(r + b, n) for x = 0, ..., n.

stat.dnbinom(s, xn, pr)
stat.pnbinom(s, xn, pr)
stat.qnbinom(p, xn, pr)

Computes the probability mass function (dnbinom), cumulative distribution function (pnbinom), and inverse cdf (qnbinom) for the negative binomial distribution with size xn >= 0 and probability of success on each trial 0 <= pr <= 1. The density is given by choose(xn + s - 1, s) * pr^xn * (1 - pr)^s for s = 0, ..., xn, where choose(xn + s - 1, s) is the binomial coefficient.

stat.dnorm(x [, mean [, sd]])
stat.pnorm(x [, mean [, sd]])
stat.qnorm(p [, mean [, sd]])

Computes the density (dnorm), cumulative distribution function (pnorm), and inverse cdf (qnorm) for the normal distribution with parameters mean and standard deviation sd >= 0. The default values for mean and sd are zero and one, respectively.

stat.dpois(s, l)
stat.ppois(s, l)
stat.qpois(p, l)

Computes the probability mass function (dpois), cumulative distribution function (ppois), and inverse cdf (qpois) for the Poisson distribution with mean l >= 0. The density is given by l^s * exp(-l) / s! for s = 0, 1, ....

stat.dt(x, df)
stat.pt(x, df)
stat.qt(p, df)

Computes the density (dt), cumulative distribution function (pt), and inverse cdf (qt) for the Student t distribution with df >= 0 degrees of freedom.

rng

Random number generation and deviates

The rng module provides random number generators based on the Mersenne-Twister generator and netlib's ranlib. Objects can be created with rng.new, copied with rng.copy, and have their seed set with rng.seed.

The following methods sample deviates from a number of probability distributions or elements in a matrix object:

      rbeta    rchisq    rexp          rf  rgamma    rnorm
    rmvnorm     runif  runifx  rdirichlet  rbinom  rnbinom
      rpois  runifint  sample     lsample

Each one of these methods are closed on a rng object and each object r has a (uservalue) table containing its methods, that is, methods for which r is an upvalue. To retrieve the method table for r, use #r. For convenience, a "class" rng object c is always created and has its methods stored in the class table so that rng.f is (#c).f.

All random deviate methods accept an optional last parameter v that must be a real matrix that serves as destination: if v is provided, v:size() deviates are generated and stored in-place at v.

Examples:

    > n = 5
    > for i = 1, n do print(i, rng.rnorm()) end
    1       -0.83608755044204
    2       0.30498104946162
    3       -1.227283848856
    4       -0.91148894025705
    5       0.28245551259444
    > r = rng.new()
    > v = (#r).rnorm(0, 1, matrix.new(n))
    > for i = 1, n do print(i, v[i]) end
    1       -0.83608755044204
    2       0.30498104946162
    3       -1.227283848856
    4       -0.91148894025705
    5       0.28245551259444

rng.new([seed])

Returns a new rng object r with seed seed. To obtain a table with methods that use r, call #r.

See also: rng, rng.seed

rng.seed([seed])

Sets the seed of rng to seed; seed can be a number or a real vector.

rng.rbeta(a, b [, dest])

Returns random deviates from the beta distribution with shape parameters a (a > 0) and b (b > 0).

See also: stat.dbeta

rng.rchisq(df [,xnonc [, dest]])

Returns random deviates from the chi-square distribution with degrees of freedom df (df > 0) and optional non-centrality parameter xnonc (xnonc >= 0, defaults to zero).

See also: stat.dchisq

rng.rexp([av [, dest]])

Returns random deviates from the exponential distribution with mean av (defaults to one).

See also: stat.dexp

rng.rf(dfn, dfd [, xnonc [, dest]])

Returns random deviates from the F distribution with degrees of freedom dfn and dfd (dfn > 0 and dfd > 0) and optional non-centrality parameter xnonc (xnonc >= 0, defaults to zero).

See also: stat.df

rng.rgamma(a [, s [, dest]])

Returns random deviates from the gamma distribution with shape a and optional rate s.

See also: stat.dgamma

rng.rnorm([m [, sd [, dest]]])

Returns random deviates from the normal distribution with mean m (defaults to zero) and standard deviation sd (sd > 0, defaults to one).

See also: stat.dnorm, rng.rmvnorm

rng.rmvnorm(m, S [, dest])

Returns random deviates from the multivariate normal distribution with mean m, a real vector of size n, and "standard deviation" S -- either a square matrix of order n, taken as the lower triangular Cholesky decomposition of the variance-covariance matrix, or a positive real vector of size n representing the square root of the diagonal of the variance-covariance matrix.

Returns S * u + m, where u is sampled from N(0, I_n).

See also: matrix.chol, rng.rnorm

rng.runif([low [, high [, dest]]])

Returns random deviates from the uniform distribution between low (defaults to zero) and high (defaults to one).

See also: stat.dunif

rng.runifx([low [, high [, dest]]])

Returns random deviates from the uniform distribution between low (defaults to zero) and high (defaults to one) with 53-bit resolution.

See also: stat.dunif

rng.rdirichlet(alpha [, dest])

Returns random deviates from the Dirichlet distribution with parameter alpha. alpha should be a real vector with positive entries.

rng.rbinom(n, p [, dest])

Returns random deviates from the binomial distribution with size n (n >= 0) and probability p (0 <= p <= 1).

See also: stat.dbinom

rng.rnbinom(n, p [, dest])

Returns random deviates from the negative binomial distribution with size n (n >= 0) and probability p (0 <= p <= 1).

See also: stat.dnbinom

rng.rpois(l [, dest])

Returns random deviates from the Poisson distribution with mean l >= 0.

See also: stat.dpois

rng.runifint([low [, high [, dest]]])

Returns random deviates from the discrete uniform distribution between low and high.

rng.sample(m [, normalized])

Returns a sample from 1 to m:size"*" with weights specified by the entries in the matrix m. If normalized is true, the entries of m are considered to be normalized, that is, it is assumed that m:sum() == 1.

rng.lsample(m [, normalized])

Returns a sample from 1 to m:size"*" with log-weights specified by the entries in the matrix m. If normalized is true, the entries of m are considered to be normalized, that is, it is assumed that m:exp():sum() == 1.

mathx

Extended math library

This module contains C99 math functions, Bessel functions (from Netlib's AMOS library), and other utilitary functions. C99 math functions are listed below:

Please refer to the man pages for details about each of these functions.

mathx.airya (z, deriv, scaling)

Returns the complex Airy function Ai(z) or its derivative if deriv is true, and a completion status. If scaling is true, returns Ai(z) * exp(zeta) where zeta = 2 / 3 * z * sqrt(z) to remove the exponential decay in -pi / 3 < arg(z) < pi / 3 and the exponential growth in pi / 3 < abs(arg(z)) < pi.

If status is true, the return is normal, that is, the computation completed without errors, underflows, or overflows; otherwise, status is an error message. If no computation is performed, nil is returned.

See also: mathx.airyb

mathx.airyb (z, deriv, scaling)

Returns the complex Airy function Bi(z) or its derivative if deriv is true, and a completion status. If scaling is true, returns Bi(z) * exp(-abs(real(zeta))) where zeta = 2 / 3 * z * sqrt(z), to remove the exponential behavior in both the left and right half planes.

If status is true, the return is normal, that is, the computation completed without errors, underflows, or overflows; otherwise, status is an error message. If no computation is performed, nil is returned.

See also: mathx.airya

mathx.besselh(nu, z, secondkind, scaling [, n])

Returns a n member sequence of complex Hankel (Bessel) functions CY(j) = H(m, nu + j - 1, z) where m = 2 if secondkind is true or m = 1 otherwise, with nonnegative orders nu + j - 1 for j = 1, ..., n and complex z ~= 0 in the cut plane -pi < arg(z) < pi. The size of the sequence n defaults to 1. A completion status is also returned.

If scaling is true, returns CY(j) * exp(-mm * z * I), where mm = 3 - 2 * m and I is the imaginary unit. The scaling removes the exponential behavior in both the upper and lower half planes.

If status is true, the return is normal, that is, the computation completed without errors, underflows, or overflows; otherwise, status is an error message. If no computation is performed, nil is returned.

See also: mathx.besseli, mathx.besselj, mathx.besselk, mathx.bessely.

mathx.besseli(nu, z, scaling [, n])

Returns a n member sequence of complex modified Bessel functions of first kind CY(j) = I(nu + j - 1, z) with nonnegative orders nu + j - 1 for j = 1, ..., n and complex z in the cut plane -pi < arg(z) < pi. The size of the sequence n defaults to 1. A completion status is also returned.

If scaling is true, returns CY(j) * exp(-abs(real(z))) to remove the exponential growth in both the left and right half planes.

If status is true, the return is normal, that is, the computation completed without errors, underflows, or overflows; otherwise, status is an error message. If no computation is performed, nil is returned.

See also: mathx.besselh, mathx.besselj, mathx.besselk, mathx.bessely.

mathx.besselj(nu, z, scaling [, n])

Returns a n member sequence of complex Bessel functions of first kind CY(j) = J(nu + j - 1, z) with nonnegative orders nu + j - 1 for j = 1, ..., n and complex z in the cut plane -pi < arg(z) < pi. The size of the sequence n defaults to 1. A completion status is also returned.

If scaling is true, returns CY(j) * exp(-abs(imag(z))) to remove the exponential growth in both the lower and upper half planes.

If status is true, the return is normal, that is, the computation completed without errors, underflows, or overflows; otherwise, status is an error message. If no computation is performed, nil is returned.

See also: mathx.besselh, mathx.besseli, mathx.besselk, mathx.bessely.

mathx.besselk(nu, z, scaling [, n])

Returns a n member sequence of complex modified Bessel functions of second kind CY(j) = K(nu + j - 1, z) with nonnegative orders nu + j - 1 for j = 1, ..., n and complex z ~= 0 in the cut plane -pi < arg(z) < pi. The size of the sequence n defaults to 1. A completion status is also returned.

If scaling is true, returns CY(j) * exp(z) to remove the exponential behavior in both the left and right half planes.

If status is true, the return is normal, that is, the computation completed without errors, underflows, or overflows; otherwise, status is an error message. If no computation is performed, nil is returned.

See also: mathx.besselh, mathx.besseli, mathx.besselj, mathx.bessely.

mathx.bessely(nu, z, scaling [, n])

Returns a n member sequence of complex Bessel functions of second kind CY(j) = Y(nu + j - 1, z) with nonnegative orders nu + j - 1 for j = 1, ..., n and complex z in the cut plane -pi < arg(z) < pi. The size of the sequence n defaults to 1. A completion status is also returned.

If scaling is true, returns CY(j) * exp(-abs(imag(z))) to remove the exponential growth in both the lower and upper half planes.

If status is true, the return is normal, that is, the computation completed without errors, underflows, or overflows; otherwise, status is an error message. If no computation is performed, nil is returned.

See also: mathx.besselh, mathx.besseli, mathx.besselj, mathx.besselk.

mathx.feq(x1, x2 [, epsilon])

Approximate equal up to precision epsilon, based on Knuth's suggestion in: Knuth, D. E. (1998), "The Art of Computer Programming. Vol. 2: Seminumerical Algorithms", Section 4.2.2.

Note that if x1 ~== x2 when feq(x1, x2) is true, then similarly

mathx.log1pe(x)

Computes log(1 + exp(x)) by avoiding calls to log and exp when abs(x) < eps.

mathx.lse(x1, x2)

Computes the "soft max" of x1 and x2, the log sum of exps log(exp(x1) + exp(x2)).

mathx.beta(x, y)
mathx.lbeta(x, y)

Computes the beta function B with arguments x and y and its logarithm (beta and lbeta respectively.) The beta function is defined by B(x, y) = integral_0^1 t^(x - 1) (1 - t)^(y - 1) dt, or, equivalently, by B(x, y) = G(x) * G(y) / G(x + y) where G is the gamma function.

See also: mathx.gamma, mathx.lgamma

mathx.digamma(x)

Computes the digamma function, the logarithmic derivative of the gamma function, at x.

See also: mathx.gamma

mathx.choose(n, k)
mathx.lchoose(n, k)

Compute binomial coefficients and their logarithms (choose and lchoose respectively`.)

License

Copyright (c) 2011 Luis Carvalho

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.