The file http://www.cs.berkeley.edu/~fateman/generic/quad-maxima.lisp
has an implementation in maxima bigfloats of Gaussian quadrature that
operates in the bfloat precision of your system. You must provide it with a
function of a single bigfloat argument (which should produce a bigfloat or
any numeric quantity that can be converted to a bigfloat), and a number of
points to use. The program will integrate the function from -1 to 1.
Choosing more points will usually get you a more accurate answer, as will
raising the fpprec value, if you have enough points.
If you want to integrate over a different range, you can do so by modifying
the function. ("change of variables")
You are invited to devise some convergence criterion for numerical
integration. There are references to papers with some suggestions.
As a part of the process, gaussianquad use a lisp program $ab_and_wts to
compute two tables: of abscissa values and weights. These can be
precomputed and (if they have sufficient precision) re-used. If this
information is available it makes integration considerably faster.
You are also invited to write a program that caches (or "memoizes") these
tables and therefore runs much faster, for a second run. These tables are
independent of the function being integrated. An example of how to do this
is given in the (lisp) version, referenced in the program.
The "main" program is written in the Macsyma language to make it more
accessible to the non-Lisp programmer. Right now that is embedded in a
Comment.
It's about 190 lines of code and comments.
Suggestions, comments, bugfixes etc are welcome.
Perhaps the lisp code (and a lisp version of the Gaussian integration should
be put in the float.lisp code.
RJF