This note explains the use of multi-precision quadrature using MPFR / lisp
on windows, Allegro CL 8.0 (possibly the free "Express" version), and
generic arithmetic. The programs could probably be converted with modest
effort to any other common lisp with a foreign function interface, though I
suspect that Allegro is particularly adept at doing this efficiently.
First, make a directory somewhere. call it "generic". Download all the
files that are mentioned below. For the case of fasl files, they
should probably be produced locally to correspond to the exact version
of Lisp that you are using. To
produce it: in Lisp, assuming the directory is c:\generic...
:cd c:/generic
:cf ga.lisp
:cf mpfr.lisp
:cf quad-fast.lisp
Next, here's what to do for running quadrature. Load the program and set the
package to mpfr
:cd c:/generic
:ld mysysquad.lisp
:pa :mpfr
(int_gs_l5bf #'exp 10) ;; use gaussian 10-point integration for -1 to 1 of
exp(x). 53 bits
(set-prec 600) ;; set precision to 600 bits
(int_gs_l5bf #'exp 10) ;; try again
(int_gs_l5bf #'exp 80) ;; try again
(time (int_gs_l5bf #'exp 80)) ;; l5 is the 5th legendre version I wrote...
This should be exp(1)-exp(-1), approx.
The program "memoizes" the weights etc to high precision, so you can use the
integration function with the same points for other integrands.
I do not include programs that attempt to achieve a particular bound by
using two or more number of nodes or precisions. The program is about 4
pages long and could be rewritten in maxima bigfloats, but would be slower.
If someone is interested in rewriting in Maxima or in lisp+bigfloats, the
source code and the documentation should be helpful. Or look up Gaussian
integration with Legendre polynomials.
The details