gsl-shell.git - gsl-shell

diff --git a/.gitignore b/.gitignore
index 997cfe37..887771ff 100644
--- a/.gitignore
+++ b/.gitignore

@@ -5,6 +5,10 @@

*.exe

*.cm[ioa]

.trash

+build

+doc/build

+www

+luajit2/src/buildvm

luajit2/lib/vmdef.lua

luajit2/src/lj_*def.h

luajit2/src/lj_vm.s

diff --git a/base.lua b/base.lua
index d0dbb5c2..73040436 100644
--- a/base.lua
+++ b/base.lua

@@ -1,8 +1,28 @@

-local cat = table.concat

-local insert = table.insert

+ -- base.lua

+ --

+ -- Copyright (C) 2009 Francesco Abbate

+ --

+ -- This program is free software; you can redistribute it and/or modify

+ -- it under the terms of the GNU General Public License as published by

+ -- the Free Software Foundation; either version 3 of the License, or (at

+ -- your option) any later version.

+ --

+ -- This program is distributed in the hope that it will be useful, but

+ -- WITHOUT ANY WARRANTY; without even the implied warranty of

+ -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

+ -- General Public License for more details.

+ --

+ -- You should have received a copy of the GNU General Public License

+ -- along with this program; if not, write to the Free Software

+ -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

+ --

-function divmod(n, p)

+local gsl = gsl or _G

+

+local cat, insert = table.concat, table.insert

+

+function gsl.divmod(n, p)

local r = n % p

return (n-r)/p, r

end

@@ -48,19 +68,18 @@ tos = function (t, maxdepth)

end

local function myprint(...)

- local n, ls = select('#', ...), {}

- for i=1, n do ls[i] = select(i, ...) end

- for i, v in ipairs(ls) do

+ local n = select('#', ...)

+ for i=1, n do

if i > 1 then io.write(', ') end

- io.write(tos(v, 3))

+ io.write(tos(select(i, ...), 3))

end

io.write('\n')

end

-echo = print

-print = myprint

+gsl.echo = print

+gsl.print = myprint

-function sequence(f, a, b)

+function gsl.sequence(f, a, b)

a, b = (b and a or 1), (b and b or a)

if not b or type(b) ~= 'number' then

error 'argument #2 should be an integer number'

@@ -76,17 +95,20 @@ end

-- take the function f and return an iterator that gives the couple (x, f(x))

-- for x going from 'xi' to 'xs' with n sampling points

-function sample(f, xi, xs, n)

+function gsl.sample(f, xi, xs, n)

local c = (xs-xi)/n

- return sequence(function(k) return xi+k*c, f(xi+k*c) end, 0, n)

+ return gsl.sequence(function(k)

+ local x = xi+k*c

+ return x, f(x)

+ end, 0, n)

end

-function ilist(f, a, b)

+function gsl.ilist(f, a, b)

local ls = {}

- for x in sequence(f, a, b) do insert(ls, x) end

+ for x in gsl.sequence(f, a, b) do insert(ls, x) end

return ls

end

-function isample(f, a, b)

- return sequence(function(i) return i, f(i) end, a, b)

+function gsl.isample(f, a, b)

+ return gsl.sequence(function(i) return i, f(i) end, a, b)

end

diff --git a/contour.lua b/contour.lua
index ccf68eff..33fbba36 100644
--- a/contour.lua
+++ b/contour.lua

@@ -17,9 +17,12 @@

-- along with this program; if not, write to the Free Software

-- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

--

+local gsl = gsl or _G

+local math = math or _G

-local insert = table.insert

-local default_color_map = color_function('redyellow', 0.9)

+local insert, sqrt, abs, divmod = table.insert, math.sqrt, math.abs, gsl.divmod

+

+local default_color_map = gsl.color_function('redyellow', 0.9)

local function reverse(ls)

local k, n = 1, #ls

@@ -29,29 +32,25 @@ local function reverse(ls)

end

end

+local function quad_approx(fl, fc, fr)

+ local a0, a1, a2 = (fr + 2*fc + fl)/4, (fr - fl)/2, (fr - 2*fc + fl)/4

+ if abs(a2) > 0.05 * abs(a1) then return end

+ local r0, r2 = (a0-a2)/a1, 2*a2/a1

+ local q = -r0 - r2*r0^2 - 2*r2^2*r0^3

+ if q >= -1 and q <= 1 then return q end

+end

+

local function root_solve(f, z0, x, y, dx0, dy0, z_eps)

local dx, dy = dx0, dy0

- local fl, fc, fr

-

- local yield = function(q) return x + q * dx, y + q * dy end

- local qeval = function(q) return f(x + q * dx, y + q * dy) - z0 end

-

- local function quad_approx()

- local a0, a1, a2 = (fr + 2*fc + fl)/4, (fr - fl)/2, (fr - 2*fc + fl)/4

- if abs(a2) > 0.05 * abs(a1) then return end

- local r0, r2 = (a0-a2)/a1, 2*a2/a1

- local q = -r0 - r2*r0^2 - 2*r2^2*r0^3

- if q >= -1 and q <= 1 then return q end

- end

for k=1,10 do

- fl, fc, fr = qeval(-1), qeval(0), qeval(1)

+ local fl, fc, fr = f(x-dx, y-dy) - z0, f(x, y) - z0, f(x+dx, y+dy) - z0

- if abs(fl) < z_eps then return yield(-1)

- elseif abs(fr) < z_eps then return yield( 1) end

+ if abs(fl) < z_eps then return x-dx, y-dy

+ elseif abs(fr) < z_eps then return x+dx, y+dy end

- local q = quad_approx()

- if q then return yield(q) else

+ local q = quad_approx(fl, fc, fr)

+ if q then return x + q * dx, y + q * dy else

dx, dy = dx/2, dy/2

local sign = (fc * fl < 0 and -1 or 1)

x, y = x + sign * dx, y + sign * dy

@@ -352,25 +351,21 @@ local function grid_create(f_, lx1, ly1, rx2, ry2, nx, ny, nlevels_or_levels, co

end

end

- local function zmap()

- for i=0, ny do

- for j=0, nx do

- local x, y = grid_point(i, j)

- local z = f(x, y)

- if not (z == z) then

- local msg = string.format('function eval at: (%g, %g) gives %g',

- x, y, z)

- error(msg)

- end

- if z < g.zmin then g.zmin = z end

- if z > g.zmax then g.zmax = z end

- g.z[index(i,j)] = z

+ for i=0, ny do

+ for j=0, nx do

+ local x, y = grid_point(i, j)

+ local z = f(x, y)

+ if not (z == z) then

+ local msg = string.format('function eval at: (%g, %g) gives %g',

+ x, y, z)

+ error(msg)

end

+ if z < g.zmin then g.zmin = z end

+ if z > g.zmax then g.zmax = z end

+ g.z[index(i,j)] = z

end

end

- zmap()

-

zlevels = {}

if type(nlevels_or_levels) == 'table' then

table.sort(nlevels_or_levels)

@@ -679,7 +674,7 @@ local function grid_create(f_, lx1, ly1, rx2, ry2, nx, ny, nlevels_or_levels, co

end

local function domain_draw_path(domid)

- local ln = path()

+ local ln = gsl.path()

local n0, n1

for _, join in ipairs(domains[domid]) do

local na, nb = curve_extrema(curves[join.id], join.direction)

@@ -716,7 +711,7 @@ local function grid_create(f_, lx1, ly1, rx2, ry2, nx, ny, nlevels_or_levels, co

end

local function curve_draw(pl, id)

- local ln = path()

+ local ln = gsl.path()

curve_add_path(ln, id, 'acw')

local tree = order_tree[id]

if tree then

@@ -745,7 +740,7 @@ local function grid_create(f_, lx1, ly1, rx2, ry2, nx, ny, nlevels_or_levels, co

end

local function grid_draw_lines(pl, col)

- local ln = path()

+ local ln = gsl.path()

for id = 1, #curves do

curve_add_path(ln, id, 'cw')

end

@@ -773,7 +768,7 @@ local function opt_gener(options, defaults)

end

end

-function contour(f, x1, y1, x2, y2, options)

+function gsl.contour(f, x1, y1, x2, y2, options)

local opt = opt_gener(options, {gridx= 40, gridy= 40, levels= 10,

colormap= default_color_map,

lines= true, show= true})

@@ -783,7 +778,7 @@ function contour(f, x1, y1, x2, y2, options)

g.find_curves()

- local p = plot()

+ local p = gsl.plot()

g.draw_regions(p)

if opt 'lines' then g.draw_lines(p, 'black') end

if opt 'show' then p:show() end

@@ -791,7 +786,7 @@ function contour(f, x1, y1, x2, y2, options)

return p

end

-function polar_contour(f, R, options)

+function gsl.polar_contour(f, R, options)

local opt = opt_gener(options, {gridx= 40, gridy= 40, levels= 10,

colormap= default_color_map,

lines= true, show= true})

@@ -801,7 +796,7 @@ function polar_contour(f, R, options)

g.find_curves()

- local p = plot()

+ local p = gsl.plot()

g.draw_regions(p)

if opt 'lines' then g.draw_lines(p, 'black') end

if opt 'show' then p:show() end

diff --git a/draw.lua b/draw.lua
index 4f1ca09c..5f84207e 100644
--- a/draw.lua
+++ b/draw.lua

@@ -1,14 +1,18 @@

-function ipath(f)

- local ln = path(f())

+local gsl = gsl or _G

+

+local floor = math and math.floor or floor

+

+function gsl.ipath(f)

+ local ln = gsl.path(f())

for x, y in f do

ln:line_to(x, y)

end

return ln

end

-function ipathp(f)

- local ln = path()

+function gsl.ipathp(f)

+ local ln = gsl.path()

local move, line = ln.move_to, ln.line_to

local function next(op)

local x, y = f()

@@ -25,34 +29,34 @@ function ipathp(f)

return ln

end

-function fxline(f, xi, xs, n)

+function gsl.fxline(f, xi, xs, n)

n = n and n or 512

- return ipath(sample(f, xi, xs, n))

+ return gsl.ipath(gsl.sample(f, xi, xs, n))

end

-function filine(f, a, b)

- return ipath(isample(f, a, b))

+function gsl.filine(f, a, b)

+ return gsl.ipath(gsl.isample(f, a, b))

end

-function xyline(x, y)

- local n = dim(x)

- local ln = path(x[1], y[1])

+function gsl.xyline(x, y)

+ local n = gsl.dim(x)

+ local ln = gsl.path(x[1], y[1])

for i=2, n do ln:line_to(x[i], y[i]) end

return ln

end

-function fxplot(f, xi, xs, color, n)

+function gsl.fxplot(f, xi, xs, color, n)

n = n and n or 512

- local p = plot()

- p:addline(ipathp(sample(f, xi, xs, n)), color)

+ local p = gsl.plot()

+ p:addline(gsl.ipathp(gsl.sample(f, xi, xs, n)), color)

p:show()

return p

end

-function fiplot(f, a, b, color)

+function gsl.fiplot(f, a, b, color)

if not b then a, b, color = 1, a, b end

- local p = plot()

- p:addline(ipathp(isample(f, a, b)), color)

+ local p = gsl.plot()

+ p:addline(gsl.ipathp(gsl.isample(f, a, b)), color)

p:show()

return p

end

@@ -65,8 +69,8 @@ local function add_square(p, lx, by, rx, ty)

p:close()

end

-function ibars(f)

- local b = path()

+function gsl.ibars(f)

+ local b = gsl.path()

local lx, ly = f()

local first = true

for rx, ry in f do

@@ -78,26 +82,22 @@ function ibars(f)

return b

end

-function segment(x1, y1, x2, y2)

- local p = path(x1, y1)

+function gsl.segment(x1, y1, x2, y2)

+ local p = gsl.path(x1, y1)

p:line_to(x2, y2)

return p

end

function rect(x1, y1, x2, y2)

- local p = path()

+ local p = gsl.path()

add_square(p, x1, y1, x2, y2)

return p

end

-function square(x0, y0, l)

- return rect(x0-l/2, y0-l/2, x0+l/2, y0+l/2)

-end

-

local bcolors = {'red', 'blue', 'green', 'magenta', 'cyan', 'yellow'}

local mcolors = {'', 'dark', 'light'}

-function rainbow(n)

+function gsl.rainbow(n)

local p = #bcolors

local q = floor((n-1)/p) % #mcolors

return mcolors[q+1] .. bcolors[(n-1) % p + 1]

@@ -109,16 +109,33 @@ local color_schema = {

darkgreen = {0.9, 0.9, 0, 0, 0.4, 0}

}

-function color_function(schema, alpha)

+function gsl.color_function(schema, alpha)

local c = color_schema[schema]

return function(a)

- return rgba(c[1] + a*(c[4]-c[1]),

- c[2] + a*(c[5]-c[2]),

- c[3] + a*(c[6]-c[3]), alpha)

+ return gsl.rgba(c[1] + a*(c[4]-c[1]),

+ c[2] + a*(c[5]-c[2]),

+ c[3] + a*(c[6]-c[3]), alpha)

end

end

-function hsl2rgb(h, s, l)

+local function HueToRgb(m1, m2, hue)

+ local v

+ hue = hue % 1

+

+ if 6 * hue < 1 then

+ v = m1 + (m2 - m1) * hue * 6

+ elseif 2 * hue < 1 then

+ v = m2

+ elseif 3 * hue < 2 then

+ v = m1 + (m2 - m1) * (2/3 - hue) * 6

+ else

+ v = m1

+ end

+

+ return v

+end

+

+function gsl.hsl2rgb(h, s, l)

local m1, m2, hue

local r, g, b

@@ -135,32 +152,15 @@ function hsl2rgb(h, s, l)

g = HueToRgb(m1, m2, h)

b = HueToRgb(m1, m2, h - 1/3)

end

- return rgb(r, g, b)

-end

-

-function HueToRgb(m1, m2, hue)

- local v

- hue = hue % 1

-

- if 6 * hue < 1 then

- v = m1 + (m2 - m1) * hue * 6

- elseif 2 * hue < 1 then

- v = m2

- elseif 3 * hue < 2 then

- v = m1 + (m2 - m1) * (2/3 - hue) * 6

- else

- v = m1

- end

-

- return v

+ return gsl.rgb(r, g, b)

end

-function hue(a)

- return hsl2rgb(a*0.7, 1, 0.6)

+function gsl.hue(a)

+ return gsl.hsl2rgb(a*0.7, 1, 0.6)

end

-function plot_lines(ln, title)

- local p = plot(title)

+function gsl.plot_lines(ln, title)

+ local p = gsl.plot(title)

for k=1, #ln do

p:addline(ln[k], rainbow(k))

end

diff --git a/gsl-shell-jit.c b/gsl-shell-jit.c
index 3dafdb59..0fe781cd 100644
--- a/gsl-shell-jit.c
+++ b/gsl-shell-jit.c

@@ -42,15 +42,18 @@ pthread_mutex_t gsl_shell_mutex[1];

static const luaL_Reg lualibs[] = {

{ "", luaopen_base },

- { LUA_LOADLIBNAME, luaopen_package },

- { LUA_TABLIBNAME, luaopen_table },

- { LUA_IOLIBNAME, luaopen_io },

- { LUA_OSLIBNAME, luaopen_os },

- { LUA_STRLIBNAME, luaopen_string },

-// { LUA_MATHLIBNAME, luaopen_math },

- { LUA_DBLIBNAME, luaopen_debug },

- { LUA_BITLIBNAME, luaopen_bit },

- { LUA_JITLIBNAME, luaopen_jit },

+ { LUA_LOADLIBNAME, luaopen_package },

+ { LUA_TABLIBNAME, luaopen_table },

+ { LUA_IOLIBNAME, luaopen_io },

+ { LUA_OSLIBNAME, luaopen_os },

+ { LUA_STRLIBNAME, luaopen_string },

+ { LUA_BITLIBNAME, luaopen_bit },

+ { LUA_JITLIBNAME, luaopen_jit },

+ { LUA_DBLIBNAME, luaopen_debug },

+#ifdef LUA_STRICT

+ { LUA_MATHLIBNAME, luaopen_math },

+ {MLUA_GSLLIBNAME, luaopen_gsl},

+#endif

{ NULL, NULL }

};

@@ -63,6 +66,7 @@ static void gsl_shell_openlibs(lua_State *L)

lua_call(L, 1, 0);

}

+#ifndef LUA_STRICT

lua_pushvalue (L, LUA_GLOBALSINDEX); /* open math in global scope */

lua_setglobal (L, LUA_MATHLIBNAME);

luaopen_math (L);

@@ -76,7 +80,7 @@ static void gsl_shell_openlibs(lua_State *L)

lua_pop (L, 1);

lua_pushnil (L); /* remove gsl */

lua_setglobal (L, MLUA_GSLLIBNAME);

-

+#endif

}

static void lstop(lua_State *L, lua_Debug *ar)

@@ -531,12 +535,7 @@ static int pmain(lua_State *L)

lua_gc(L, LUA_GCSTOP, 0); /* stop collector during initialization */

gsl_shell_openlibs(L); /* open libraries */

lua_gc(L, LUA_GCRESTART, -1);

-

- dolibrary (L, "base");

- dolibrary (L, "integ");

- dolibrary (L, "igsl");

- dolibrary (L, "draw");

-

+ dolibrary (L, "gslext");

s->status = handle_luainit(L);

if (s->status != 0) return 0;

script = collectargs(argv, &has_i, &has_v, &has_e);

diff --git a/gsl-shell.c b/gsl-shell.c
index ecb9d6cb..f3c31ba2 100644
--- a/gsl-shell.c
+++ b/gsl-shell.c

@@ -85,9 +85,9 @@ static const luaL_Reg gshlibs[] = {

{LUA_IOLIBNAME, luaopen_io},

{LUA_OSLIBNAME, luaopen_os},

{LUA_STRLIBNAME, luaopen_string},

+ {LUA_DBLIBNAME, luaopen_debug},

#ifdef LUA_STRICT

{LUA_MATHLIBNAME, luaopen_math},

- {LUA_DBLIBNAME, luaopen_debug},

{MLUA_GSLLIBNAME, luaopen_gsl},

#endif

{NULL, NULL}

@@ -493,16 +493,8 @@ static int pmain (lua_State *L) {

if (argv[0] && argv[0][0]) progname = argv[0];

lua_gc(L, LUA_GCSTOP, 0); /* stop collector during initialization */

gsl_shell_openlibs (L); /* open libraries */

-

lua_gc(L, LUA_GCRESTART, 0);

-

- dolibrary (L, "base");

- dolibrary (L, "integ");

-#ifndef LUA_STRICT

- dolibrary (L, "igsl");

- dolibrary (L, "draw");

-#endif

-

+ dolibrary (L, "gslext");

s->status = handle_luainit(L);

if (s->status != 0) return 0;

script = collectargs(argv, &has_i, &has_v, &has_e);

diff --git a/gslext.lua b/gslext.lua
new file mode 100644
index 00000000..ef77265b
--- /dev/null
+++ b/gslext.lua

@@ -0,0 +1,6 @@

+-- load initialization files for GSL Shell

+

+require('base')

+require('matrix')

+require('integ')

+require('draw')

diff --git a/integ.lua b/integ.lua
index 61f8e256..e7e97489 100644
--- a/integ.lua
+++ b/integ.lua

@@ -1,5 +1,7 @@

-- integr module follows

+local gsl = gsl or _G

+

local function get_spec(spec, key, mtype)

local defaults = {eps_abs= 1e-6, eps_rel= 1e-6, limit= 512,

rule= 'SINGULAR'}

@@ -17,7 +19,7 @@ local function type_check(v, mtype, name)

end

end

-function integ(spec)

+function gsl.integ(spec)

local limit_min = 512

local f = spec.f

@@ -59,15 +61,15 @@ function integ(spec)

if (inttype == 'awf') then

cspec.eps_abs = cspec.eps_abs / 2

cspec.a, cspec.b = 0, 0

- local lint, lerr = integ_raw(f, 'awfl', cspec, w)

- local uint, uerr = integ_raw(f, 'awfu', cspec, w)

+ local lint, lerr = gsl.integ_raw(f, 'awfl', cspec, w)

+ local uint, uerr = gsl.integ_raw(f, 'awfu', cspec, w)

return lint + uint, lerr + uerr

end

- return integ_raw(f, inttype, cspec, w)

+ return gsl.integ_raw(f, inttype, cspec, w)

end

error('cannot calculate indefinite integral with this weight')

end

- return integ_raw(f, inttype, cspec)

+ return gsl.integ_raw(f, inttype, cspec)

end

for i,v in ipairs(pts) do

@@ -91,7 +93,7 @@ function integ(spec)

else

error('unknown wight type')

end

- return integ_raw(f, inttype, cspec, w)

+ return gsl.integ_raw(f, inttype, cspec, w)

end

if not spec.adaptive then

@@ -115,5 +117,5 @@ function integ(spec)

inttype = 'ag'

end

- return integ_raw(f, inttype, cspec)

+ return gsl.integ_raw(f, inttype, cspec)

end

diff --git a/luajit2/src/Makefile b/luajit2/src/Makefile
index 9b78e3c0..2bf099f4 100644
--- a/luajit2/src/Makefile
+++ b/luajit2/src/Makefile

@@ -8,6 +8,11 @@

# Copyright (C) 2005-2010 Mike Pall. See Copyright Notice in luajit.h

##############################################################################

+GSH_BASE_DIR = ../..

+

+include $(GSH_BASE_DIR)/makeconfig

+include $(GSH_BASE_DIR)/makeflags

+

MAJVER= 2

MINVER= 0

RELVER= 0

@@ -31,7 +36,6 @@ CC= gcc

# unwinding are not affected -- the assembler part has frame unwind

# information and GCC emits it where needed (x64) or with -g (see CCDEBUG).

# CCOPT= -O2 -fomit-frame-pointer

-CCOPT= -g

# Use this if you want to generate a smaller binary (but it's slower):

#CCOPT= -Os -fomit-frame-pointer

# Note: it's no longer recommended to use -O3 with GCC 4.x.

@@ -177,7 +181,9 @@ TARGET_ASHLDFLAGS= $(LDOPTIONS) $(TARGET_XSHLDFLAGS) $(TARGET_SHLDFLAGS)

TARGET_ALIBS= $(TARGET_XLIBS) $(LIBS) $(TARGET_LIBS)

ifneq (,$(findstring stack-protector,$(shell $(TARGET_CC) -dumpspecs)))

- TARGET_XCFLAGS+= -fno-stack-protector

+ ifeq (,$(findstring no-stack-protector,$(CFLAGS)))

+ TARGET_XCFLAGS+= -fno-stack-protector

+ endif

endif

ifeq (,$(findstring __i386__,$(shell echo __i386__ | $(TARGET_CC) -P -E -)))

diff --git a/luajit2/src/luaconf.h b/luajit2/src/luaconf.h
index 4506b5fb..fde2770b 100644
--- a/luajit2/src/luaconf.h
+++ b/luajit2/src/luaconf.h

@@ -6,6 +6,9 @@

#ifndef luaconf_h

#define luaconf_h

+#define QUOTEME_(x) #x

+#define QUOTEME(x) QUOTEME_(x)

+

#include <limits.h>

#include <stddef.h>

@@ -26,30 +29,45 @@

** In Windows, any exclamation mark ('!') in the path is replaced by the

** path of the directory of the executable file of the current process.

*/

-#define LUA_LDIR "!\\lua\\"

+ #ifdef GSL_SHELL_LUA

+ #define STDLIB_NAME "gsl-shell"

+ #else

+ #define STDLIB_NAME "lua"

+ #endif

+#define LUA_LDIR "!\\" STDLIB_NAME "\\"

#define LUA_CDIR "!\\"

#define LUA_PATH_DEFAULT \

".\\?.lua;" LUA_LDIR"?.lua;" LUA_LDIR"?\\init.lua;"

#define LUA_CPATH_DEFAULT \

".\\?.dll;" LUA_CDIR"?.dll;" LUA_CDIR"loadall.dll"

#else

-#define LUA_ROOT "/usr/local/"

-#define LUA_LDIR LUA_ROOT "share/lua/5.1/"

-#define LUA_CDIR LUA_ROOT "lib/lua/5.1/"

+ #ifdef GSL_SHELL_LUA

+ #define STDLIB_NAME "gsl-shell"

+ #else

+ #define STDLIB_NAME "lua/5.1"

+ #endif

+#define LUA_ROOT_Q QUOTEME(LUA_ROOT) "/"

+#define LUA_LDIR LUA_ROOT_Q "share/" STDLIB_NAME "/"

+#define LUA_CDIR LUA_ROOT_Q "lib/" STDLIB_NAME "/"

#ifdef LUA_XROOT

#define LUA_JDIR LUA_XROOT "share/luajit-2.0.0-beta5/"

#define LUA_XPATH \

- ";" LUA_XROOT "share/lua/5.1/?.lua;" LUA_XROOT "share/lua/5.1/?/init.lua"

-#define LUA_XCPATH LUA_XROOT "lib/lua/5.1/?.so;"

+ ";" LUA_XROOT "share/" STDLIB_NAME "/?.lua;" LUA_XROOT "share/" STDLIB_NAME "/?/init.lua"

+#define LUA_XCPATH LUA_XROOT "lib/" STDLIB_NAME "/?.so;"

#else

-#define LUA_JDIR LUA_ROOT "share/luajit-2.0.0-beta5/"

+#define LUA_JDIR LUA_ROOT_Q "share/luajit-2.0.0-beta5/"

#define LUA_XPATH

#define LUA_XCPATH

#endif

-#define LUA_PATH_DEFAULT \

+ #ifdef GSL_SHELL_LUA

+ #define LUA_PATH_DEFAULT "./?.lua;" LUA_LDIR"?.lua;" LUA_LDIR"?/init.lua"

+ #define LUA_CPATH_DEFAULT "./?.so;" LUA_CDIR"?.so;" LUA_CDIR"loadall.so"

+ #else

+ #define LUA_PATH_DEFAULT \

"./?.lua;" LUA_JDIR"?.lua;" LUA_LDIR"?.lua;" LUA_LDIR"?/init.lua" LUA_XPATH

-#define LUA_CPATH_DEFAULT \

+ #define LUA_CPATH_DEFAULT \

"./?.so;" LUA_CDIR"?.so;" LUA_XCPATH LUA_CDIR"loadall.so"

+ #endif

#endif

/* Environment variable names for path overrides and initialization code. */

diff --git a/makeflags b/makeflags
index 54306887..a75be9e4 100644
--- a/makeflags
+++ b/makeflags

@@ -9,6 +9,6 @@ ifeq ($(strip $(DEBUG)), yes)

CXXFLAGS = -g -Wall

else

# on windows do not use -fomit-frame-pointer

- CFLAGS = -O3 -fomit-frame-pointer -ffast-math -Wall

- CXXFLAGS = -Os -ffast-math -fno-rtti -Wall

+ CFLAGS = -O3 -fno-stack-protector -fomit-frame-pointer -ffast-math -Wall

+ CXXFLAGS = -Os -fno-stack-protector -ffast-math -fno-rtti -Wall

endif

diff --git a/igsl.lua b/matrix.lua
index 0e0d421c..42350d51 100644
--- a/igsl.lua
+++ b/matrix.lua

@@ -1,5 +1,5 @@

- -- igsl.lua

+ -- matrix.lua

--

-- Copyright (C) 2009 Francesco Abbate

--

@@ -18,12 +18,15 @@

-- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

--

-local cat = table.concat

-local insert = table.insert

-local fmt = string.format

+local gsl = gsl or _G

+local math = math or _G

+

+local cat, insert, fmt = table.concat, table.insert, string.format

+

+local sqrt, abs = math.sqrt, math.abs

local function matrix_f_set(m, f)

- local r, c = dim(m)

+ local r, c = gsl.dim(m)

for i = 1, r do

for j = 1, c do

local z = f(i, j)

@@ -33,8 +36,8 @@ local function matrix_f_set(m, f)

return m

end

-function matrix_reduce(m, f, accu)

- local r, c = dim(m)

+function gsl.matrix_reduce(m, f, accu)

+ local r, c = gsl.dim(m)

for i = 1, r do

for j = 1, c do

accu = f(accu, m:get(i, j))

@@ -44,20 +47,20 @@ function matrix_reduce(m, f, accu)

end

local function tostring_eps(z, eps)

- local a, b = real(z), imag(z)

- local f = |x| fmt('%g', x)

- local s = abs(a) > eps and f(a) or ''

+ local a, b = gsl.real(z), gsl.imag(z)

+ local s = abs(a) > eps and fmt('%g', a) or ''

if b > eps then

local sign = (s == '' and '' or '+')

- s = s .. fmt('%s%si', sign, f(b))

+ s = s .. fmt('%s%gi', sign, b)

elseif b < -eps then

- s = s .. fmt('-%si', f(-b))

+ s = s .. fmt('-%gi', -b)

end

return s == '' and '0' or s

end

local function matrix_from_table(ctor, t)

- return matrix_f_set(ctor(#t, #t[1]), |i,j| t[i][j])

+ local r, c = #t, #t[1]

+ return matrix_f_set(ctor(r, c), function(i,j) return t[i][j] end)

end

local function vector_from_table(ctor, t)

@@ -66,99 +69,125 @@ local function vector_from_table(ctor, t)

return v

end

-function vector(t)

+function gsl.vector(t)

return vector_from_table(new, t)

end

-function cvector(t)

+function gsl.cvector(t)

return vector_from_table(cnew, t)

end

-function matrix(t)

+function gsl.matrix(t)

return matrix_from_table(new, t)

end

-function cmatrix(t)

+function gsl.cmatrix(t)

return matrix_from_table(cnew, t)

end

-function matrix_to_string(m)

+local function padstr(s, w)

+ return fmt('%s%s', string.rep(' ', w - #s), s)

+end

+

+local function matrix_to_string(m)

local eps = m:norm() * 1e-8

local fwidth = function(w, val)

local ln = # tostring_eps(val, eps)

return (ln > w and ln or w)

end

- local width = matrix_reduce(m, fwidth, 0)

- local pad = |s| string.rep(' ', width - #s) .. s

- local r, c = dim(m)

+ local width = gsl.matrix_reduce(m, fwidth, 0)

+ local r, c = gsl.dim(m)

local lines = {}

for i=1,r do

local ln = {}

for j=1,c do

- insert(ln, pad(tostring_eps(m:get(i,j), eps)))

+ insert(ln, padstr(tostring_eps(m:get(i,j), eps), width))

end

- insert(lines, '[ ' .. cat(ln, ' ') .. ' ]')

+ insert(lines, fmt('[ %s ]', cat(ln, ' ')))

end

return cat(lines, '\n')

end

-function tr(m)

- local r, c = dim(m)

- return new(c, r, |i,j| m:get(j,i))

+local function csqr(z)

+ local r, i = gsl.real(z), gsl.imag(z)

+ return r*r + i*i

+end

+

+function gsl.tr(m)

+ local r, c = gsl.dim(m)

+ return gsl.new(c, r, function(i,j) return m:get(j,i) end)

end

-function hc(m)

- local r, c = dim(m)

- return cnew(c, r, |i,j| conj(m:get(j,i)))

+function gsl.hc(m)

+ local r, c = gsl.dim(m)

+ return gsl.cnew(c, r, function(i,j) return gsl.conj(m:get(j,i)) end)

end

-function diag(v)

- local n = dim(v)

- return new(n, n, |i,j| i == j and v:get(i,1) or 0)

+function gsl.diag(v)

+ local n = gsl.dim(v)

+ return gsl.new(n, n, function(i,j) return i == j and v:get(i,1) or 0 end)

end

-function unit(n)

- return new(n, n, |i,j| i == j and 1 or 0)

+function gsl.unit(n)

+ return gsl.new(n, n, function(i,j) return i == j and 1 or 0 end)

end

-function matrix_norm(m)

- local sq = matrix_reduce(m, |p, z| p + real(z)*real(z) + imag(z)*imag(z), 0)

- return sqrt(sq)

+local function matrix_norm(m)

+ local r, c = gsl.dim(m)

+ local s = 0

+ for i=1, r do

+ for j=1, c do

+ s = s + csqr(m:get(i,j))

+ end

+ end

+ return sqrt(s)

end

-function matrix_column (m, c)

- local r = dim(m)

+local function matrix_column (m, c)

+ local r = gsl.dim(m)

return m:slice(1, c, r, 1)

end

-function matrix_row (m, r)

- local _, c = dim(m)

+local function matrix_row (m, r)

+ local _, c = gsl.dim(m)

return m:slice(r, 1, 1, c)

end

-function matrix_rows(m)

- local r, c = dim(m)

- return sequence(|i| m:slice(i, 1, 1, c), r)

+local function matrix_rows(m)

+ local r, c = gsl.dim(m)

+ return gsl.sequence(function(i) m:slice(i, 1, 1, c) end, r)

end

-function set(d, s)

- matrix_f_set(d, |i,j| s:get(i,j))

+function gsl.set(d, s)

+ local r, c = gsl.dim(m)

+ local dset, sget = d.set, s.get

+ for i=1, r do

+ for j=1, c do

+ dset(i, j, sget(i, j))

+ end

+ end

end

-function null(m)

- matrix_f_set(m, |i,j| 0)

+function gsl.null(m)

+ local r, c = gsl.dim(m)

+ local mset = m.set

+ for i=1, r do

+ for j=1, c do

+ mset(i, j, 0)

+ end

+ end

end

-function fset(m, f)

+function gsl.fset(m, f)

matrix_f_set(m, f)

end

local function add_matrix_method(s, m)

- Matrix[s] = m

- cMatrix[s] = m

+ gsl.Matrix[s] = m

+ gsl.cMatrix[s] = m

end

-function ode_iter(s, t0, y0, t1, h, tstep)

+function gsl.ode_iter(s, t0, y0, t1, h, tstep)

s:set(t0, y0, h)

return function()

local t, y = s.t, s.y

@@ -176,7 +205,7 @@ local function hc_reduce(hc, f, accu)

end

local function hc_print(hc)

- local eps = 1e-8 * hc_reduce(hc, |p,z| p + z*conj(z), 0)

+ local eps = 1e-8 * hc_reduce(hc, function(p,z) return p + csqr(z) end, 0)

local f = function(p, z)

insert(p, fmt('%6i: %s', #p, tostring_eps(z, eps)))

return p

@@ -184,19 +213,19 @@ local function hc_print(hc)

return cat(hc_reduce(hc, f, {}), '\n')

end

-FFT_hc_mixed_radix.__tostring = hc_print

-FFT_hc_radix2.__tostring = hc_print

+gsl.FFT_hc_mixed_radix.__tostring = hc_print

+gsl.FFT_hc_radix2.__tostring = hc_print

-ODE.iter = ode_iter

-cODE.iter = ode_iter

+gsl.ODE.iter = ode_iter

+gsl.cODE.iter = ode_iter

local function add_matrix_meta_method(key, method)

local m, mt

- m = new(1,1)

+ m = gsl.new(1,1)

mt = getmetatable(m)

mt[key] = method

- m = cnew(1,1)

+ m = gsl.cnew(1,1)

mt = getmetatable(m)

mt[key] = method

end

@@ -208,9 +237,9 @@ add_matrix_method('col', matrix_column)

add_matrix_method('row', matrix_row)

add_matrix_method('rows', matrix_rows)

-function linmodel(f, x)

- local p, n = #f(x[1]), dim(x)

- local A = new(n, p)

+function gsl.linmodel(f, x)

+ local p, n = #f(x[1]), gsl.dim(x)

+ local A = gsl.new(n, p)

for k=1,n do

local y = f(x[k])

for j=1,p do A:set(k, j, y[j]) end

@@ -218,18 +247,18 @@ function linmodel(f, x)

return A

end

-function linfit(gener, x, y, w)

- local X = linmodel(gener, x)

- local c, cov = mlinear(X, y, w)

+function gsl.linfit(gener, x, y, w)

+ local X = gsl.linmodel(gener, x)

+ local c, cov = gsl.mlinear(X, y, w)

local f = function(xe)

- local xs = vector(gener(xe))

- return prod(xs, c)[1]

+ local xs = gsl.vector(gener(xe))

+ return gsl.prod(xs, c)[1]

end

return f, c

end

local function nlinfitwrapf(fmodel, x, y)

- local n = dim(y)

+ local n = gsl.dim(y)

return function(p, f, J)

for k=1, n do

local ym = fmodel(p, x[k], J and J:row(k))

@@ -239,10 +268,10 @@ local function nlinfitwrapf(fmodel, x, y)

end

-function nlinfit(f, x, y, p0)

- local N = dim(y)

- local P = dim(p0)

- local s = nlfsolver {fdf= nlinfitwrapf(f, x, y), n= N, p0= p0}

+function gsl.nlinfit(f, x, y, p0)

+ local N = gsl.dim(y)

+ local P = gsl.dim(p0)

+ local s = gsl.nlfsolver {fdf= nlinfitwrapf(f, x, y), n= N, p0= p0}

s:run()

- return |x| f(s.p, x), s.p

+ return function(x) return f(s.p, x) end, s.p

end