[ANN] Numeric Lua 0.1

[Luis Carvalho <carvalho@...>]

Ola',

I'd like to announce the availability of Numeric Lua. Numeric Lua is a
numerical (duh!) package for Lua. It includes support for complex numbers,
multidimensional matrices, random number generation and special functions.
Most of the routines are simple wrappers for stable and well known libraries
from Netlib (http://www.netlib.org).

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

Some features/highlights:

- Matrices have a "type" attribute -- general, symmetric, triangular, or
  positive definite -- for faster (and/or more suitable) operations. A
  standard interface to BLAS/LAPACK routines allows the use of optimized
  routines from ATLAS or MKL, for example.

- Matrices, complex numbers, and RNG's (random number generators) can be seen
  as objects. Matrices have some functional facilities for faster, C-sided
  computations.

- RNG's use the fast Mersenne Twister generator as an "engine" (uniform
  deviates) and Netlib's Ranlib for other deviates (normal, Poisson, gamma,
  ...)

- Special functions include Airy, Bessel, gamma related (gamma, lgamma,
  psigamma, factorial, beta, lbeta, ...), exponentials (err function, expm1,
  log1p, ...), CDF's, PDF's and others.

Drawbacks:

- Only tested and known to run on Lua 5.1 work6.

- No documentation at all. :(

- Very, very beta. Be warned! :)


Numeric Lua can be found at:

    http://luaforge.net/projects/numlua

Finally, a quick test drive:

$ lua
Lua 5.1 (work6)  Copyright (C) 1994-2005 Tecgraf, PUC-Rio
> -- let's generate some random numbers (fast!)
> require 'rng'
> r = rng(os.time())
> for i=1,5 do print(r:unif()) end -- numbers can vary, of course :)
0.47017181466799
0.63861816644203
0.55608703091275
0.68107319541741
0.24100433581043
>
> -- now, let's add some math to spice up
> require 'spfun'
> = spfun.log1p(1e-16) -- keep precision
1e-16
> = math.log(1+1e-16)
0
> = spfun.isinf(1/0)
true
> = spfun.isnan(0/0)
true
> = spfun.qnorm(0.025, 0, 1) -- normal quantiles, mean = 0, var = 1
-1.9599639845401
> x = 0.025; print(spfun.pnorm(spfun.qnorm(x, 0, 1), 0, 1) == x)
true
>
> -- complex numbers
> require 'complex'
> = complex(1,2) + complex.i
1 + 3i
> x = complex.sqrt(-2)
> = complex.exp(x) -- cos(sqrt(2)) + sin(sqrt(2))*i
0.15594369476537 + 0.98776594599274i
>
> -- matrices
> require 'matrix'
> x = matrix(2,3,4) -- 2 x 3 x 4 matrix
> for k, i, j in x:entries() do x[k][i][j] = k * (i + j) end
> matrix.foreach(x, print) -- list x
...
> -- let's try the type system
> x = matrix.zeros(1000) -- 1000 x 1000 matrix
> x:apply(function(i, j) return i + j end) -- x[i][j] = i+j
> = x:gettype() -- x is general
G
> -- let's profile transpose(x) * x
> local t = os.clock(); print((*x) * x); print(os.clock() - t)
matrix: 0xb586001c
1.32
> x:settype() -- guess x's type
> = x:gettype() -- x is symmetric
S
> local t = os.clock(); print((*x) * x); print(os.clock() - t)
1.28
> -- a bit faster... let's compute the inner product directly
> local t = os.clock(); print(x:inner()); print(os.clock() - t)
0.66

Cheers,
luis.

-- 
A mathematician is a device for turning coffee into theorems.
        -- P. Erdos 

-- 
Luis Carvalho
Applied Math PhD Student - Brown University
PGP Key: E820854A <carvalho@dam.brown.edu>

Attachment: pgpBl5dYqTYUT.pgp
Description: PGP signature