Skip to content

Latest commit

 

History

History
126 lines (96 loc) · 3.94 KB

gsl-ffi.org

File metadata and controls

126 lines (96 loc) · 3.94 KB

Using emacs-ffi with the Gnu Scientific Library

First you define the library. We use the GNU Scientific library here (https://www.gnu.org/software/gsl/doc/html/index.html)

(require 'ffi)

(define-ffi-library gsl "libgsl")

Simple functions

We use a Bessel function here.

(define-ffi-function gsl-sf-bessel-J0 "gsl_sf_bessel_J0"
  (:double "Function value") 		;return value
  ((:double x "x")) ; args
  gsl
  "Regular cylindrical Bessel function of zeroth order, J_0(x)")
(gsl-sf-bessel-J0 5.0)

The macro gives a reasonable docstring.

(describe-function 'gsl-sf-bessel-J0)

Integration function

This is a more complex example where we integrate a function. This block sets up some GSL integration functions, and a struct.

(define-ffi-function gsl-integration-workspace-alloc "gsl_integration_workspace_alloc"
  (:pointer "gsl_integration_workspace")
  ((:size_t n))
  gsl
  "This function allocates a workspace sufficient to hold n
double precision intervals, their integration results and error
estimates. One workspace may be used multiple times as all
necessary reinitialization is performed automatically by the
integration routines.")

(define-ffi-function gsl-integration-workspace-free "gsl_integration_workspace_free"
  (:void)
  ((:pointer *w "gsl_integration_workspace"))
  gsl
  "This function frees the memory associated with the workspace w.")

(define-ffi-function gsl-integration-qags "gsl_integration_qags"
  (:int "Integral result")
  ((:pointer *f "gsl_function")
   (:double a) (:double b)
   (:double epsabs) (:double epsrel)
   (:size_t limit) (:pointer *w "gsl-integration-workspace")
   (:pointer *result "double") (:pointer *abserr "double"))
  gsl
  "This function applies the Gauss-Kronrod 21-point integration
rule adaptively until an estimate of the integral of f over (a,b)
is achieved within the desired absolute and relative error
limits, epsabs and epsrel. The results are extrapolated using the
epsilon-algorithm, which accelerates the convergence of the
integral in the presence of discontinuities and integrable
singularities. The function returns the final approximation from
the extrapolation, result, and an estimate of the absolute error,
abserr. The subintervals and their results are stored in the
memory provided by workspace. The maximum number of subintervals
is given by limit, which may not exceed the allocated size of the
workspace.")

(define-ffi-struct gsl_function (function :type :pointer) (params :type :pointer))

Here is a working example. We use it to integrate the function $f(x) = log x / sqrt(x)$ from 0 to 1. This is a hard one since it starts at log 0 which is undefined, and 1 / 0, which is also undefined. But, the integral converges to -4.

There is a lot going on in here. We integrate a lisp function (that is turned into a pointer with ffi-make-closure. The result is stored in a float we have to allocate, and then retrieve with ffi--mem-ref.

(defun ff (x params) (/ (log x) (sqrt x)))

(let* ((*w (gsl-integration-workspace-alloc 1000))
       (F (ffi--define-struct gsl_function))
       (cif (ffi--prep-cif :double [:double :pointer]))
       (*f (ffi-make-closure cif #'ff))
       (a 0.0)
       (b 1.0)
       (epsabs 0.0)
       (epsrel 1e-7)
       (limit 1000)
       (result (ffi-allocate :double))
       (abserr (ffi-allocate :double)))
  (setf (gsl_function-function F) *f)
  (gsl-integration-qags F a b epsabs epsrel limit *w result abserr)
  (ffi--mem-ref result :double))

Success!