Gaussian RNG?, numerics in general...



>>>>> "Richard" == Richard Fateman <fateman@cs.berkeley.edu> writes:

    Richard> Liam M. Healy wrote:


    >> 
    >> This is a great idea, but sadly, FFIs are not standard between CL
    >> implementations. 

I should add: As Erik Naggum once observed, the hard part isn't foreign _function_
interface, it's foreign (usually meaning C) _data_structure_
interface, and this has been my experience.  I would also add that
there's a foreign signal interface problem as well (see below).

    Richard> I think Franz Inc has a program that reads header files and generates
    Richard> the FFI necessary  (at some level of effectiveness < 100%).

cbind.  It requires a lot of hand work, but it's better than nothing.
Last I checked, it only runs on specific platforms not including
Linux.  I have seen similar things mentioned
  http://www.bricoworks.com/~moore/cparse/index.html
  http://www.zipcon.net/~briand/ffigen.html
but I don't know anything about them.  

    >> 
    >> I would very much like to see GSL (http://sources.redhat.com/gsl/) be
    >> accessible from CL; it seems to have just about everything one could
    >> want in a sci/math library, and is actively being developed. 


    Richard> This is probably a negative from the point of view of support.  That is,
    Richard> active development ...

Huh?  There seems to be quite a bit of progress on it.  I have posted
questions to their mailing list and received help, so there is some
informal support.  Jim Amundson recommends GSL too, so that counts for
something here :-).


    >> This
    >> library might well be a useful resource for Maxima to build on.  I
    >> have used this from ACL on a function-by-function basis but the hard
    >> part is building the interface. 
    >> 

    Richard> Generally yes, but a sufficiently "functional" view helps.
    Richard> GMP does mess up its parameters somewhat, and it does some of
    Richard> its own memory management. but I think it is really OK.


    Richard> Nonstandard FFI is a problem;  I am still trying to get the Franz
    Richard> people to do something that makes sense for maxima and for them so
    Richard> that we don't need so many lisps.  This may be impossible because
    Richard> some people are committed to their own lisps.  And their own front
    Richard> ends.
    Richard>    At least we don't have pitched battles (so far).

    Richard> I would be interested in any material you might have for ACL and linking
    Richard> to GSL.
    Richard> I'm still learning about the "new" foreign function stuff.

OK.  Here's some of my code, to access the elliptic functions:

#+allegro
(ignore-errors (load "/usr/local/lib/libgslcblas.so"))

#+allegro
(ignore-errors (load "/usr/local/lib/libgsl.so"))

#+allegro
(foreign-functions:def-foreign-call
    ;;; See http://sources.redhat.com/gsl/ref/gsl-ref_7.html#SEC94
    (jacobian-elliptic-function-int "gsl_sf_elljac_e")
    ((u :double)
     (m :double)
     (sn (:array :double)) (cn (:array :double)) (dn (:array :double)))
  :convention :stdcall
  :returning :int)

(defvar sn (make-array 1 :element-type 'double-float))
(defvar cn (make-array 1 :element-type 'double-float))
(defvar dn (make-array 1 :element-type 'double-float))

(defun jacobian-elliptic-cn (u m)
  "The Jacobian elliptic function cn(u,m)."
  (jacobian-elliptic-function-int u m sn cn dn)
  (aref cn 0))

(defun jacobian-elliptic-sn (u m)
  "The Jacobian elliptic function sn(u,m)."
  (jacobian-elliptic-function-int u m sn cn dn)
  (aref sn 0))

Now here there are a couple examples of the problems we face in
linking with a foreign library.  First, notice that I have a more
Lispy (functional) layer on top of the FFI, because I don't want to
have to cobble up the double float locations for return values on each
call.  Unless the auto generators are really fancy, I think you have
to do this by hand.  Second, if m>1, the elliptic functions are
undefined.  GSL signals an error, but it is handled by default by
termination, so ACL crashes.  GSL is set up so that you can define
your own error handler (in C, of course), and I ask Franz how to set
up a foreign signal handling that CL would handle intelligently (at
the very least, not crashing).  They gave me some instruction but I
never implemented it; it sounded a bit iffy.

Having said that, I think using the existing library is the only way
to go.  The suggestions to reimplement in CL, while having a certain
appeal, are just not practical.  We need to leverage the existing
resources, and reimpmelementaion would take a huge amount of time and
effort, and on an ongoing basis as well.

I also think we're getting away from maxima-specific things and ought
to bring this to comp.lang.lisp or one of the other general CL
discussions.  There ought to be broader interest than just Maxima.


-- 
Liam Healy
Liam.Healy@nrl.navy.mil