Gaussian RNG?, numerics in general...



Sigh. I had really hoped to postpone this conversation until I felt we
had fixed a few more underpinnings of Maxima (like the packages!), but
the topic has suddenly heated up. Before I go on to agree with Richard
and to make a concrete proposal, let me say this:

Some (if not much, or even all) of the numerical code in Maxima is in
really bad shape. Only after Ray Toy started fixing up the bessel code
did I notice that the Maxima test suite has been dutifully verifying
that 
	jn(3,4) = 0.1320342093678868 ;
	j0(1) = 0.7651977429270344 ;
	bessel(2,3) = 0.1289432631185581-5.717744444379852e-08*%i;
even though only the first six digits or so of each expression is
correct. Yes, Maxima always gets the same answer to 16 or so digits, but
most of the digits are wrong. bessel(2,3) should be real, but we've been
verifying that it gets the wrong imaginary part to sixteen digits!

Writing reliable numerical software is a hard job. I think it is
completely unreasonable to think we can go about improving the numerical
part of Maxima one function at a time. Ray's work is fine as a stopgap
measure, but we need to think about a real long-term solution.

On Mon, 2002-02-04 at 11:42, Richard Fateman wrote:
> If someone has a program in C or Fortran which takes as input
> a few integers and/or floats and returns, after some substantial
> computation,  a float, an integers, or an array of these, then it should
> be possible to link this to maxima easily via foreign function
> calls, and not ever worry about lisp's speed.

Yes, absolutely. I'm worried about maintenance even more than speed. If
bugs are identified and fixed in the original C/Fortran/whatever code,
we want to be able to take advantage of them immediately.

For me, the biggest issue is the foreign function interfaces in the
various Lisps on various platforms. Unfortunately, I don't know anything
about it! I hope someone on the list can enlighten me as to some of the
issues.

> For a random number generator you might even produce an array
> of a whole bunch of these from one call.
> 
> An idea which really capitalizes on the fact that we are
>   using lisp and maxima as open source is that we can patch in
> any old stuff from a "subset domain".  Thus any numerical
> code should be simple.  Octave, a Matlab clone, could be
> swallowed whole, I think.

Interesting. Actually, I had a different idea for combining Maxima and
Octave, but I think we differ primarily in details. Maybe we can have
that conversation later...

> For example, over the weekend I tried out (successfully) loading
> gmp 3.0  (gnu multiple precision) into a lisp. Interestingly, I
> did NOT need the source for gmp. I loaded a dll.  All i needed
> was the documentation for the entry points.
> 
> 
> Other features which come free from other peoples' lisp libraries
> include network socket code, html generation.
> 
> The major hassle in my view in doing this is that other peoples' code
> may not work on all relevant platforms, or may be subject to random
> fluctuations.  E.g. incompatible changes in interfaces, legal status,
> etc.

These are all things to worry about. I propose we use the GNU Scientific
Library (GSL) to replace {X} in Maxima, where the set {X} includes at
least special functions and random number generators. GSL has the
following features:
	- For most of what it does, it is the best such code I know of.
	- It is actively maintained.
	- It is documented.
	- It has been ported to a large number of platforms.
	- Its license is completely compatible with ours.
A number of smart people have done a *large* amount of work on GSL. We
should take advantage of it.

--Jim