diff --git a/calculus/numerical-integration-with-tables.lisp b/calculus/numerical-integration-with-tables.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..75584db8a84b4b0186715617b2ffab9b84abf5d1
--- /dev/null
+++ b/calculus/numerical-integration-with-tables.lisp
@@ -0,0 +1,254 @@
+;; Numerical integration techniques that require tables
+;; Liam Healy 2009-04-04 15:24:05EDT 
+;; Time-stamp: <2009-04-04 18:09:47EDT numerical-integration-with-tables.lisp>
+;; $Id: $
+
+(in-package :gsl)
+
+;;; /usr/include/gsl/gsl_integration.h
+
+;;;;****************************************************************************
+;;;; QAWS adaptive integration for singular functions
+;;;;****************************************************************************
+
+(defmobject qaws-table
+    "gsl_integration_qaws_table"
+  ((alpha :double) (beta :double) (mu :int) (nu :int))
+  "table for QAWS numerical integration method"
+  :documentation			; FDL
+  "Make and initialize a table for the QAWS
+   adaptive integration method for singular functions.
+   It a singular weight function W(x) with the parameters (alpha, beta, mu, nu),
+          W(x) = (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x)
+   where alpha > -1, beta > -1, and mu = 0, 1, nu = 0, 1. The
+   weight function can take four different forms depending on the
+   values of mu and nu,
+          W(x) = (x-a)^alpha (b-x)^beta                   (mu = 0, nu = 0)
+          W(x) = (x-a)^alpha (b-x)^beta log(x-a)          (mu = 1, nu = 0)
+          W(x) = (x-a)^alpha (b-x)^beta log(b-x)          (mu = 0, nu = 1)
+          W(x) = (x-a)^alpha (b-x)^beta log(x-a) log(b-x) (mu = 1, nu = 1)
+   The singular points (a,b) do not have to be specified until the
+   integral is computed, where they are the endpoints of the integration
+   range."
+  :initialize-suffix "set"
+  ;; This defines the reinitialize-instance but causes the make-qaws-table
+  ;; to initialize twice.
+  :initialize-args ((alpha :double) (beta :double) (mu :int) (nu :int)))
+
+(defmfun integration-QAWS
+    (function a b alpha beta mu nu
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (table (make-qaws-table alpha beta mu nu))
+	      (limit 1000) (workspace (make-integration-workspace limit)))
+  "gsl_integration_qaws"
+  ((callback :pointer)
+   (a :double) (b :double)
+   ((mpointer table) :pointer)
+   (absolute-error :double) (relative-error :double)
+   (limit sizet) ((mpointer workspace) :pointer)
+   (result :double) (abserr :double))
+  :callbacks
+  (callback gsl-function nil (function :double (:input :double) :slug))
+  :callback-dynamic (nil (function))
+  :documentation			; FDL
+  "Compute the integral of the function f(x) over the interval (a,b)
+  with the singular weight function (x-a)^alpha (b-x)^beta
+  log^mu (x-a) log^nu (b-x). The parameters of the weight
+  function (alpha, beta, mu, nu) are used to make the default table. The
+  integral is
+    I = int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x).
+  The adaptive bisection algorithm of QAG is used. When a
+  subinterval contains one of the endpoints then a special 25-point
+  modified Clenshaw-Curtis rule is used to control the
+  singularities. For subintervals which do not include the endpoints
+  an ordinary 15-point Gauss-Kronrod integration rule is used.")
+
+;;;;****************************************************************************
+;;;; QAWO adaptive integration for oscillatory functions
+;;;;****************************************************************************
+
+(cffi:defcenum integrate-sine-cosine
+  :cosine :sine)
+
+(defmobject qawo-table
+    "gsl_integration_qawo_table"
+  ((omega :double) (L :double) (trig integrate-sine-cosine) (n sizet))
+  "table for QAWO numerical integration method"
+  :documentation			; FDL
+  "Make a table describing a sine or cosine weight function W(x) with
+  the parameters (omega, L),
+          W(x) = sin(omega x)
+          W(x) = cos(omega x)
+  The parameter L must be the length of the interval over which the
+  function will be integrated L = b - a. The choice of sine or cosine
+  is made with the parameter trig which should be one of :cosine or :sine.
+  This makes a table of the trigonometric coefficients required in the
+  integration process. The parameter n determines the number of levels
+  of coefficients that are computed. Each level corresponds to one
+  bisection of the interval L, so that n levels are sufficient for
+  subintervals down to the length L/2^n.  An error of class
+  'table-limit-exceeded is signalled if the number of levels is
+  insufficient for the requested accuracy."
+  :initialize-suffix "set"
+  ;; This defines the reinitialize-instance but causes the make-qawo-table
+  ;; to initialize twice.
+  :initialize-args ((omega :double) (L :double) (trig integrate-sine-cosine)))
+
+(defmfun integration-QAWO
+    (function a omega L trig n
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (table (make-qawo-table omega L trig n))
+	      (limit 1000) (workspace (make-integration-workspace limit)))
+  "gsl_integration_qawo"
+  ((callback :pointer)
+   (a :double)
+   (absolute-error :double) (relative-error :double)
+   (limit sizet) ((mpointer workspace) :pointer)
+   ((mpointer table) :pointer)
+   (result :double) (abserr :double))
+  :callbacks
+  (callback gsl-function nil (function :double (:input :double) :slug))
+  :callback-dynamic (nil (function))
+  :documentation			; FDL
+  "Use an adaptive algorithm to compute the integral of f over (a,b)
+   with the weight function \sin(\omega x) or \cos(\omega x) defined
+   by the table wf,
+          I = \int_a^b dx f(x) sin(omega x)
+          I = \int_a^b dx f(x) cos(omega x)
+
+   The results are extrapolated using the epsilon-algorithm to
+   accelerate the convergence of the integral. 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.
+
+   Those subintervals with large widths d where d\omega > 4 are
+   computed using a 25-point Clenshaw-Curtis integration rule, which
+   handles the oscillatory behavior. Subintervals with a small
+   widths where d\omega < 4 are computed using a 15-point
+   Gauss-Kronrod integration.")
+
+;;;;****************************************************************************
+;;;; QAWF adaptive integration for Fourier integrals
+;;;;****************************************************************************
+
+(defmfun integration-QAWF
+    (function a omega L trig n
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (table (make-qawo-table omega L trig n))
+	      (limit 1000)
+	      (workspace (make-integration-workspace limit))
+	      (cycle-workspace (make-integration-workspace limit)))
+  "gsl_integration_qawf"
+  ((callback :pointer)
+   (a :double)
+   (absolute-error :double)
+   (limit sizet)
+   ((mpointer workspace) :pointer) ((mpointer cycle-workspace) :pointer)
+   ((mpointer table) :pointer)
+   (result :double) (abserr :double))
+  :callbacks
+  (callback gsl-function nil (function :double (:input :double) :slug))
+  :callback-dynamic (nil (function))
+  :documentation			; FDL
+  "This function attempts to compute a Fourier integral of the
+   function f over the semi-infinite interval [a,+\infty).
+
+              I = \int_a^{+\infty} dx f(x) sin(omega x)
+              I = \int_a^{+\infty} dx f(x) cos(omega x)
+
+    The parameter \omega and choice of \sin or \cos is taken from the
+    table wf (the length L can take any value, since it is overridden
+    by this function to a value appropriate for the fourier
+    integration). The integral is computed using the QAWO algorithm
+    over each of the subintervals,
+
+              C_1 = [a, a + c]
+              C_2 = [a + c, a + 2 c]
+              ... = ...
+              C_k = [a + (k-1) c, a + k c]
+
+    where c = (2 floor(|\omega|) + 1) \pi/|\omega|. The width c is
+    chosen to cover an odd number of periods so that the contributions
+    from the intervals alternate in sign and are monotonically
+    decreasing when f is positive and monotonically decreasing. The
+    sum of this sequence of contributions is accelerated using the
+    epsilon-algorithm.
+
+    This function works to an overall absolute tolerance of
+    abserr. The following strategy is used: on each interval C_k the
+    algorithm tries to achieve the tolerance
+
+              TOL_k = u_k abserr
+
+    where u_k = (1 - p)p^{k-1} and p = 9/10. The sum of the geometric
+    series of contributions from each interval gives an overall
+    tolerance of abserr.
+
+    If the integration of a subinterval leads to difficulties then the
+    accuracy requirement for subsequent intervals is relaxed,
+
+              TOL_k = u_k max(abserr, max_{i<k}{E_i})
+
+    where E_k is the estimated error on the interval C_k.
+
+    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. The integration over each subinterval uses the memory
+    provided by cycle_workspace as workspace for the QAWO algorithm.")
+
+;;;;****************************************************************************
+;;;; Examples and unit test
+;;;;****************************************************************************
+
+(defun integration-test-f456 (x)
+  (if (zerop x)
+      0.0d0
+      ;; Note this function is different than the C case; it takes the
+      ;; absolute value of the argument and then transfers the sign of
+      ;; x to the answer.
+      (* (signum x) (log (abs x)))))
+
+(defun integration-test-f457 (x)
+  (if (zerop x)
+      0.0d0
+      (/ (sqrt x))))
+
+(defun integration-test-f458 (x)
+  (if (zerop x)
+      0.0d0
+      (expt (1+ (expt (log x) 2)) -2)))
+
+;; Tests from gsl-1.11/integration/test.c
+;; Functions defined in gsl-1.11/integration/tests.c
+(save-test numerical-integration
+ (integration-QAWS
+  'integration-test-f458 0.0d0 1.0d0 0.0d0 0.0d0 1 0 0.0d0 1.0d-7)
+ (integration-QAWS
+  'integration-test-f458 0.0d0 1.0d0 -0.5d0 -0.3d0 0 0 0.0d0 1.0d-7)
+ (integration-QAWS
+  'integration-test-f458 0.0d0 1.0d0 -0.5d0 -0.3d0 1 0 0.0d0 1.0d-7)
+ (integration-QAWS
+  'integration-test-f458 0.0d0 1.0d0 -0.5d0 -0.3d0 0 1 0.0d0 1.0d-7)
+ (integration-QAWS
+  'integration-test-f458 0.0d0 1.0d0 -0.5d0 -0.3d0 1 1 0.0d0 1.0d-7)
+ ;; Thought test.c has n=1000, I find that this causes a floating
+ ;; point overflow in making the table.  Setting to 100 works and
+ ;; gives the result given in test.c.
+ (integration-QAWO
+  'integration-test-f456 0.0d0 (* 10 pi) 1.0d0 :sine 100 0.0d0 1.0d-7)
+ ;; For the negative interval, it was necessary to have the function
+ ;; takethe absolute value of the argument and then transfer the sign of
+ ;;  x to the answer.
+ (integration-QAWO
+  'integration-test-f456 0.0d0 (* 10 pi) -1.0d0 :sine 100 0.0d0 1.0d-7)
+ (integration-QAWF
+  'integration-test-f457 0.0d0 (/ pi 2) 1.0d0 :cosine 1000 1.0d-7))
diff --git a/calculus/numerical-integration.lisp b/calculus/numerical-integration.lisp
index 7740919ead9e48daef9cddf955a8ffc7b474dbf0..fcb9371279ab0771503d015aa42e3b12c2775d3f 100644
--- a/calculus/numerical-integration.lisp
+++ b/calculus/numerical-integration.lisp
@@ -1,20 +1,32 @@
 ;; Numerical integration
 ;; Liam Healy, Wed Jul  5 2006 - 23:14
-;; Time-stamp: <2009-03-31 22:09:43EDT numerical-integration.lisp>
+;; Time-stamp: <2009-04-04 19:35:50EDT numerical-integration.lisp>
 ;; $Id$
 
-;;; To do: QAWS, QAWO, QAWF, more tests
-
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_integration.h
 
+;;;;****************************************************************************
+;;;; Default error values
+;;;;****************************************************************************
+
+(export '(*default-absolute-error* *default-relative-error*))
+
+(defparameter *default-absolute-error* 1.0d-5
+  "The default absolute error used in numerical integration.")
+
+(defparameter *default-relative-error* 0.0d0
+  "The default relative error used in numerical integration.")
+
 ;;;;****************************************************************************
 ;;;; QNG non-adaptive Gauss-Kronrod integration
 ;;;;****************************************************************************
 
 (defmfun integration-QNG
-    (function a b &optional (absolute-error 1.0d0) (relative-error 1.0d0))
+    (function a b
+	      &optional (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*))
   ;; Set absolute-error and relative-error to 1 because it apparently doesn't matter
   ;; what these are if they are too large, it will do a minimum number
   ;; of points anyway.
@@ -26,7 +38,7 @@
   :callbacks
   (callback gsl-function nil (function :double (:input :double) :slug))
   :callback-dynamic (nil (function))
-  :documentation				       ; FDL
+  :documentation			; FDL
   "Apply the Gauss-Kronrod 10-point, 21-point, 43-point and
    87-point integration rules in succession until an estimate of the
    integral of f over (a,b) is achieved within the desired
@@ -53,9 +65,11 @@
   :gauss41 :gauss51 :gauss61)
 
 (defmfun integration-QAG
-    (function a b method limit 
-	      &optional (absolute-error 1.0d0) (relative-error 1.0d0)
-	      (workspace (make-integration-workspace limit)))
+    (function a b method
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (limit 1000) (workspace (make-integration-workspace limit)))
   ;; Set absolute-error and relative-error to 1 because it apparently doesn't matter
   ;; what these are if they are too large, it will do a minimum number
   ;; of points anyway.
@@ -92,9 +106,10 @@
 ;;;;****************************************************************************
 
 (defmfun integration-QAGS
-    (function a b limit
-	      &optional (absolute-error 1.0d0) (relative-error 1.0d0)
-	      (workspace (make-integration-workspace limit)))
+    (function a b 
+	      &optional (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (limit 1000) (workspace (make-integration-workspace limit)))
   "gsl_integration_qags"
   ((callback :pointer)
    (a :double) (b :double)
@@ -122,12 +137,14 @@
 ;;;;****************************************************************************
 
 (defmfun integration-QAGP
-    (function points limit
-	      &optional (absolute-error 1.0d0) (relative-error 1.0d0)
-	      (workspace (make-integration-workspace limit)))
+    (function points
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (limit 1000) (workspace (make-integration-workspace limit)))
   "gsl_integration_qagp"
   ((callback :pointer)
-   ((mpointer points) :pointer) ((dim0 points) sizet)
+   ((c-pointer points) :pointer) ((dim0 points) sizet)
    (absolute-error :double) (relative-error :double) (limit sizet)
    ((mpointer workspace) :pointer) (result :double) (abserr :double))
   :inputs (points)
@@ -151,9 +168,11 @@
 ;;;;****************************************************************************
 
 (defmfun integration-QAGi
-    (function limit
-	      &optional (absolute-error 1.0d0) (relative-error 1.0d0)
-	      (workspace (make-integration-workspace limit)))
+    (function 
+     &optional
+     (absolute-error *default-absolute-error*)
+     (relative-error *default-relative-error*)
+     (limit 1000) (workspace (make-integration-workspace limit)))
   "gsl_integration_qagi"
   ((callback :pointer)
    (absolute-error :double) (relative-error :double) (limit sizet)
@@ -173,9 +192,11 @@
    this case a lower-order rule is more efficient.")
 
 (defmfun integration-QAGiu
-    (function a limit
-	      &optional (absolute-error 1.0d0) (relative-error 1.0d0)
-	      (workspace (make-integration-workspace limit)))
+    (function a
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (limit 1000) (workspace (make-integration-workspace limit)))
   "gsl_integration_qagiu"
   ((callback :pointer) (a :double)
    (absolute-error :double) (relative-error :double) (limit sizet)
@@ -191,9 +212,11 @@
    and then integrated using the QAGS algorithm.")
 
 (defmfun integration-QAGil
-    (function b limit
-	      &optional (absolute-error 1.0d0) (relative-error 1.0d0)
-	      (workspace (make-integration-workspace limit)))
+    (function b
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (limit 1000) (workspace (make-integration-workspace limit)))
   "gsl_integration_qagil"
   ((callback :pointer) (b :double)
    (absolute-error :double) (relative-error :double) (limit sizet)
@@ -213,9 +236,11 @@
 ;;;;****************************************************************************
 
 (defmfun integration-QAWC
-    (function a b c limit
-	      &optional (absolute-error 1.0d0) (relative-error 1.0d0)
-	      (workspace (make-integration-workspace limit)))
+    (function a b c 
+	      &optional
+	      (absolute-error *default-absolute-error*)
+	      (relative-error *default-relative-error*)
+	      (limit 1000) (workspace (make-integration-workspace limit)))
   "gsl_integration_qawc"
   ((callback :pointer)
    (a :double) (b :double) (c :double)
@@ -241,7 +266,97 @@
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
+;;; CCL 1.2 returns a complex number when the exponent is a double
+;;; float even if it's positive whole number value.  This is a
+;;; workaround.
+#+ccl
+(defun nn-expt (base power)
+  (if (= power (floor power))
+      (expt base (floor power))
+      (expt base power)))
+
+(defun integration-test-f1 (alpha)
+  (lambda (x) (* (expt x alpha) (log (/ x)))))
+
+(defun integration-test-f3 (alpha)
+  (lambda (x) (cos (* (expt 2 alpha) (sin x)))))
+
+(defun integration-test-f11 (alpha)
+  (lambda (x) (#+ccl nn-expt #-ccl expt (log (/ x)) (1- alpha))))
+
+(defun integration-test-f15 (alpha)
+  (lambda (x) (* x x (exp (* (- (expt 2 (- alpha))) x)))))
+
+(defun integration-test-f16 (alpha)
+  (lambda (x)
+    (cond
+      ((and (= alpha 1.0d0) (zerop x)) 1.0d0)
+      ((and (> alpha 1.0d0) (zerop x)) 0.0d0)
+      (t (/ (#+ccl nn-expt #-ccl expt x (1- alpha))
+	    (expt (1+ (* 10 x)) 2))))))
+
+(defun integration-test-f454 (x)
+  (* (expt x 3) (* (log (abs (* (- (expt x 2) 1.0d0) (- (expt x 2) 2.0d0)))))))
+
+(defun integration-test-f455 (x)
+  (/ (log x) (1+ (* 100 x x))))
+
+(defun integration-test-f459 (x)
+  (/ (+ (* 5.0d0 (expt x 3)) 6.0d0)))
+
+(defun integration-test-myfn1 (x)
+  (exp (- (- x) (expt x 2))))
+
+(defun integration-test-myfn2 (alpha)
+  (lambda (x) (exp (* alpha x))))
+
 (save-test numerical-integration
-  (integration-qng 'sin 0.0d0 pi)
-  (integration-QAG 'sin 0.0d0 pi :gauss15 20)
-  (integration-QAG 'sin 0.0d0 pi :gauss21 40))
+ (integration-qng 'sin 0.0d0 pi)
+ (integration-QAG 'sin 0.0d0 pi :gauss15 20)
+ (integration-QAG 'sin 0.0d0 pi :gauss21 40)
+ ;; Tests from gsl-1.11/integration/test.c
+ ;; Functions defined in gsl-1.11/integration/tests.c
+ (integration-QNG (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.1d0 0.0d0)
+ (integration-QNG (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.1d0 0.0d0)
+ (integration-QNG (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.0d0 1.0d-9)
+ (integration-QNG (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.0d0 1.0d-9)
+ (integration-QNG (integration-test-f3 1.3d0) 0.3d0 2.71d0 0.0d0 1d-12)
+ (integration-QNG (integration-test-f3 1.3d0) 2.71d0 0.3d0 0.0d0 1d-12)
+ (integration-QNG (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.0d0 1.0d-13)
+ (integration-QNG (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.0d0 1.0d-13)
+ (integration-QNG (integration-test-f1 -0.9d0) 0.0d0 1.0d0 0.0d0 1.0d-3) ; error
+ (integration-QNG (integration-test-f1 -0.9d0) 1.0d0 0.0d0 0.0d0 1.0d-3) ; error
+ (integration-QAG
+  (integration-test-f1 2.6d0) 0.0d0 1.0d0 :gauss15 0.0d0 1.0d-10 1000)
+ (integration-QAG
+  (integration-test-f1 2.6d0) 1.0d0 0.0d0 :gauss15 0.0d0 1.0d-10 1000)
+ (integration-QAG
+  (integration-test-f1 2.6d0) 0.0d0 1.0d0 :gauss21 1.0d-14 0.0d0 1000)
+ (integration-QAG
+  (integration-test-f1 2.6d0) 1.0d0 0.0d0 :gauss21 1.0d-14 0.0d0 1000)
+ (integration-QAG			; roundoff error
+  (integration-test-f3 1.3d0) 0.3d0 2.71d0 :gauss31 1.0d-14 0.0d0 1000)
+ (integration-QAG			; roundoff error
+  (integration-test-f3 1.3d0) 2.71d0 0.3d0 :gauss31 1.0d-14 0.0d0 1000)
+ (integration-QAG			; singularity error
+  (integration-test-f16 2.0d0) -1.0d0 1.0d0 :gauss51 1.0d-14 0.0d0 1000)
+ (integration-QAG			; singularity error
+  (integration-test-f16 2.0d0) 1.0d0 -1.0d0 :gauss51 1.0d-14 0.0d0 1000)
+ (integration-QAG			; iteration limit error
+  (integration-test-f16 2.0d0) -1.0d0 1.0d0 :gauss61 1.0d-14 0.0d0 3)
+ (integration-QAG			; iteration limit error
+  (integration-test-f16 2.0d0) 1.0d0 -1.0d0 :gauss61 1.0d-14 0.0d0 3)
+ (integration-QAGS (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.0d0 1d-10 1000)
+ (integration-QAGS (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.0d0 1d-10 1000)
+ (integration-QAGS (integration-test-f11 2.0d0) 1.0d0 1000.0d0 1d-7 0.0d0 1000)
+ (integration-QAGS (integration-test-f11 2.0d0) 1000.0d0 1.0d0 1d-7 0.0d0 1000)
+ (integration-QAGiu 'integration-test-f455 0.0d0 0.0d0 1d-3 1000)
+ (integration-QAGiu (integration-test-f15 5.0d0) 0.0d0 0.0d0 1d-7 1000)
+ (integration-QAGiu (integration-test-f16 1.0d0) 99.9d0 1d-7 0.0d0 1000)
+ (integration-QAGi 'integration-test-myfn1 1.0d-7 0.0d0 1000)
+ (integration-QAGil (integration-test-myfn2 1.0d0) 1.0d0 1.0d-7 0.0d0 1000)
+ (integration-QAGp
+  'integration-test-f454
+  (copy (vector 0.0d0 1.0d0 (sqrt 2.0d0) 3.0d0) 'vector-double-float)
+  0.0d0 1.0d-3 1000)
+ (integration-QAWc 'integration-test-f459 -1.0d0 5.0d0 0.0d0 0.0d0 1.0d-3 1000))
diff --git a/documentation/index.html b/documentation/index.html
index a3b59b56a4e98f4d8f79c1f034763e8ad0816e4f..cacee1448a41be43b48634bce517a378855dce26 100644
--- a/documentation/index.html
+++ b/documentation/index.html
@@ -96,7 +96,7 @@ complex conjugate scalar product</a> of two complex vectors of length 3:
       #2m(49.27d0 -13.49d0 32.5d0 42.73d0 -17.24d0 43.31d0))
 #C(-2940.2118d0 1861.9380999999998d0)
 </pre>
-<p>There are over 1200 examples available from within GSLL with the
+<p>There are over 1500 examples available from within GSLL with the
 function <code>examples</code>.  These examples also serve
 as a test suite for GSLL.
 
@@ -334,6 +334,7 @@ histogram histogram2d histogram-pdf histogram2d-pdf
 basis-spline chebyshev hankel wavelet wavelet-workspace
 random-number-generator quasi-random-number-generator discrete-random
 polynomial-complex-workspace integration-workspace
+qaws-table qawo-table
 eigen-symm eigen-symmv eigen-herm eigen-hermv
 eigen-nonsymm eigen-nonsymmv eigen-gensymm eigen-gensymmv
 eigen-gen eigen-genv
@@ -385,7 +386,7 @@ and arrays used internally or for function return.
 <!-- Created: Feb 25 2005 -->
 <!-- hhmts start -->
 <small>
-Time-stamp: <2009-04-04 09:58:19EDT index.html>
+Time-stamp: <2009-04-04 18:01:52EDT index.html>
 </small>
 <!-- hhmts end -->
  </div>
diff --git a/gsll.asd b/gsll.asd
index da8150f298a9a68e3aae8fd167d312d308c701a8..fe69a055156a5063d90105118955f0c33068c393 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2009-04-04 09:55:40EDT gsll.asd>
+;; Time-stamp: <2009-04-04 17:56:17EDT gsll.asd>
 ;; $Id$
 
 (asdf:defsystem "gsll"
@@ -178,6 +178,8 @@
 	    :depends-on (init data random)
 	    :components
 	    ((:file "numerical-integration")
+	     (:file "numerical-integration-with-tables"
+		    :depends-on (numerical-integration))
 	     (:file "monte-carlo")
 	     (:file "numerical-differentiation")))
    (:module ordinary-differential-equations
diff --git a/init/defmfun-single.lisp b/init/defmfun-single.lisp
index 971d47d0a5babaf7d948ea6771ce338205ef7cea..9492ddc5cd918e7f9c3452d9e4f9c8c2e7c85330 100644
--- a/init/defmfun-single.lisp
+++ b/init/defmfun-single.lisp
@@ -1,6 +1,6 @@
 ;; Helpers that define a single GSL function interface
 ;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
-;; Time-stamp: <2009-04-01 21:36:13EDT defmfun-single.lisp>
+;; Time-stamp: <2009-04-04 18:47:47EDT defmfun-single.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -99,12 +99,12 @@
 		 (append before after
 			 (callback-symbol-set
 			  callback-dynamic callbacks (first callback-dynamic-variables))
-			 (let ((auxstart (position '&aux arglist)))
-			   ;; &aux bindings are checked
+			   ;; &optional/&key/&aux defaults are checked
+			 (let ((auxstart (after-llk arglist)))
 			   (when auxstart
 			     (apply
 			      'append
-			      (mapcar 'rest (subseq arglist (1+ auxstart))))))))))))
+			      (mapcar 'rest (remove-if 'atom auxstart)))))))))))
 	     (first callback-dynamic-variables))
 	   ,@(when documentation (list documentation))
 	   ,(funcall body-maker name arglist gsl-name c-arguments key-args))
diff --git a/init/forms.lisp b/init/forms.lisp
index e257ddfada6fb16bf0ec20f13075cb602c7ec981..0b2b4916d8c521de76258d975e9478e0c0890407 100644
--- a/init/forms.lisp
+++ b/init/forms.lisp
@@ -1,6 +1,6 @@
 ;; Lisp forms
 ;; Liam Healy 2009-03-07 15:49:25EST forms.lisp
-;; Time-stamp: <2009-04-02 22:34:03EDT forms.lisp>
+;; Time-stamp: <2009-04-04 19:08:43EDT forms.lisp>
 
 (in-package :gsl)
 
@@ -40,3 +40,20 @@
       (when pos
 	(nth pos cats)))))
 
+(defun after-llk (arglist)
+  "The portion of the arglist from the first llk on."
+  (loop for elt in arglist
+       with seen = nil
+       when (or seen (member elt *defmfun-llk*))
+       do (setf seen t)
+       when seen collect elt))
+
+#|
+;;; Oddly, this gives a warning in SBCL but works the same.
+(defun after-llk (arglist)
+  "The portion of the arglist from the first llk on."
+  (when (intersection *defmfun-llk* arglist)
+    (subseq arglist
+	    (some (lambda (itm) (position itm arglist)) *defmfun-llk*))))
+|#
+
diff --git a/tests/numerical-integration.lisp b/tests/numerical-integration.lisp
index 09047741b34959185213f39371d2ab598a81f6f9..ab0c2b1cd75b423ed04881590d47236d71fb63df 100644
--- a/tests/numerical-integration.lisp
+++ b/tests/numerical-integration.lisp
@@ -2,17 +2,220 @@
 
 (in-package :gsl)
 
-(LISP-UNIT:DEFINE-TEST NUMERICAL-INTEGRATION
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 2.0d0 2.220446049250313d-14 21)
-                        (MULTIPLE-VALUE-LIST
-                         (INTEGRATION-QNG 'SIN 0.0d0 PI)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 2.0d0 2.220446049250313d-14)
-                        (MULTIPLE-VALUE-LIST
-                         (INTEGRATION-QAG 'SIN 0.0d0 PI :GAUSS15 20)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 2.0d0 2.220446049250313d-14)
-                        (MULTIPLE-VALUE-LIST
-                         (INTEGRATION-QAG 'SIN 0.0d0 PI :GAUSS21 40))))
+(LISP-UNIT:DEFINE-TEST
+    NUMERICAL-INTEGRATION
+  (let ((lisp-unit:*epsilon* (* 5.0d7 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST -0.18927518534894017d0 1.1291337120150095d-8)
+     (MULTIPLE-VALUE-LIST
+      (INTEGRATION-QAWS 'INTEGRATION-TEST-F458 0.0d0 1.0d0
+			0.0d0 0.0d0 1 0 0.0d0 1.d-7))))
+  (let ((lisp-unit:*epsilon* (* 5.0d7 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST 0.9896686656601706d0 5.888032496373315d-8)
+     (MULTIPLE-VALUE-LIST
+      (INTEGRATION-QAWS 'INTEGRATION-TEST-F458 0.0d0 1.0d0
+			-0.5d0 -0.3d0 0 0 0.0d0 1.d-7))))
+  (let ((lisp-unit:*epsilon* (* 5.0d7 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST -0.36366794705865396d0 2.851348775255562d-8)
+     (MULTIPLE-VALUE-LIST
+      (INTEGRATION-QAWS 'INTEGRATION-TEST-F458 0.0d0 1.0d0
+			-0.5d0 -0.3d0 1 0 0.0d0 1.d-7))))
+  (let ((lisp-unit:*epsilon* (* 5.0d7 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST -1.9114892533634094d0 9.854016804653183d-9)
+     (MULTIPLE-VALUE-LIST
+      (INTEGRATION-QAWS 'INTEGRATION-TEST-F458 0.0d0 1.0d0
+			-0.5d0 -0.3d0 0 1 0.0d0 1.d-7))))
+  (let ((lisp-unit:*epsilon* (* 5.0d7 double-float-epsilon)))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST 0.31599228628110465d0 2.33618347181736d-8)
+     (MULTIPLE-VALUE-LIST
+      (INTEGRATION-QAWS 'INTEGRATION-TEST-F458 0.0d0 1.0d0
+			-0.5d0 -0.3d0 1 1 0.0d0 1.d-7))))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.12813684839916742d0 6.875028324415666d-12)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAWO 'INTEGRATION-TEST-F456 0.0d0
+		      (* 10 PI) 1.0d0 :SINE 100 0.0d0
+		      1.d-7)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.12813684839916742d0 6.875028324415666d-12)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAWO 'INTEGRATION-TEST-F456 0.0d0
+		      (* 10 PI) -1.0d0 :SINE 100 0.0d0
+		      1.d-7)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9999999999279803d0 1.556289974669056d-8)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAWF 'INTEGRATION-TEST-F457 0.0d0
+		      (/ PI 2) 1.0d0 :COSINE 1000 1.d-7)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 2.0d0 2.220446049250313d-14 21)
+   (MULTIPLE-VALUE-LIST (INTEGRATION-QNG 'SIN 0.0d0 PI)))
+  (LISP-UNIT:ASSERT-ERROR 'TYPE-ERROR
+			  (INTEGRATION-QAG 'SIN 0.0d0 PI
+					   :GAUSS15 20))
+  (LISP-UNIT:ASSERT-ERROR 'TYPE-ERROR
+			  (INTEGRATION-QAG 'SIN 0.0d0 PI
+					   :GAUSS21 40))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07716049379303085d0 9.424302194248481d-8 21)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F1 2.6d0) 0.0d0
+		     1.0d0 0.1d0 0.0d0)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.07716049379303085d0 9.424302194248481d-8 21)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F1 2.6d0) 1.0d0
+		     0.0d0 0.1d0 0.0d0)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07716049382706505d0 2.666891413688592d-12 43)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F1 2.6d0) 0.0d0
+		     1.0d0 0.0d0 1.d-9)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.07716049382706505d0 2.666891413688592d-12 43)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F1 2.6d0) 1.0d0
+		     0.0d0 0.0d0 1.d-9)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.7238969575482962d0 1.2776768895200564d-14 43)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F3 1.3d0) 0.3d0
+		     2.71d0 0.0d0 1.d-12)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.7238969575482962d0 1.2776768895200564d-14 43)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F3 1.3d0) 2.71d0
+		     0.3d0 0.0d0 1.d-12)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.0771604938271603d0 8.566535680046933d-16 87)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F1 2.6d0) 0.0d0
+		     1.0d0 0.0d0 1.d-13)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.0771604938271603d0 8.566535680046933d-16 87)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QNG (INTEGRATION-TEST-F1 2.6d0) 1.0d0
+		     0.0d0 0.0d0 1.d-13)))
+  (LISP-UNIT:ASSERT-ERROR 'FAILURE-TO-REACH-TOLERANCE
+			  (INTEGRATION-QNG
+			   (INTEGRATION-TEST-F1 -0.9d0)
+			   0.0d0 1.0d0 0.0d0 0.001d0))
+  (LISP-UNIT:ASSERT-ERROR 'FAILURE-TO-REACH-TOLERANCE
+			  (INTEGRATION-QNG
+			   (INTEGRATION-TEST-F1 -0.9d0)
+			   1.0d0 0.0d0 0.0d0 0.001d0))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07716049382715855d0 6.679384889070553d-12)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAG (INTEGRATION-TEST-F1 2.6d0) 0.0d0
+		     1.0d0 :GAUSS15 0.0d0 1.d-10 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.07716049382715855d0 6.679384889070553d-12)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAG (INTEGRATION-TEST-F1 2.6d0) 1.0d0
+		     0.0d0 :GAUSS15 0.0d0 1.d-10 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07716049382716048d0 5.181082604031411d-15)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAG (INTEGRATION-TEST-F1 2.6d0) 0.0d0
+		     1.0d0 :GAUSS21 1.d-14 0.0d0 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.07716049382716048d0 5.181082604031411d-15)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAG (INTEGRATION-TEST-F1 2.6d0) 1.0d0
+		     0.0d0 :GAUSS21 1.d-14 0.0d0 1000)))
+  (LISP-UNIT:ASSERT-ERROR 'ROUNDOFF-FAILURE
+			  (INTEGRATION-QAG
+			   (INTEGRATION-TEST-F3 1.3d0)
+			   0.3d0 2.71d0 :GAUSS31 1.d-14
+			   0.0d0 1000))
+  (LISP-UNIT:ASSERT-ERROR 'ROUNDOFF-FAILURE
+			  (INTEGRATION-QAG
+			   (INTEGRATION-TEST-F3 1.3d0)
+			   2.71d0 0.3d0 :GAUSS31 1.d-14
+			   0.0d0 1000))
+  (LISP-UNIT:ASSERT-ERROR 'SINGULARITY
+			  (INTEGRATION-QAG
+			   (INTEGRATION-TEST-F16 2.0d0)
+			   -1.0d0 1.0d0 :GAUSS51 1.d-14
+			   0.0d0 1000))
+  (LISP-UNIT:ASSERT-ERROR 'SINGULARITY
+			  (INTEGRATION-QAG
+			   (INTEGRATION-TEST-F16 2.0d0)
+			   1.0d0 -1.0d0 :GAUSS51 1.d-14
+			   0.0d0 1000))
+  (LISP-UNIT:ASSERT-ERROR 'EXCEEDED-MAXIMUM-ITERATIONS
+			  (INTEGRATION-QAG
+			   (INTEGRATION-TEST-F16 2.0d0)
+			   -1.0d0 1.0d0 :GAUSS61 1.d-14
+			   0.0d0 3))
+  (LISP-UNIT:ASSERT-ERROR 'EXCEEDED-MAXIMUM-ITERATIONS
+			  (INTEGRATION-QAG
+			   (INTEGRATION-TEST-F16 2.0d0)
+			   1.0d0 -1.0d0 :GAUSS61 1.d-14
+			   0.0d0 3))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07716049382715791d0 2.216394969317974d-12)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGS (INTEGRATION-TEST-F1 2.6d0) 0.0d0
+		      1.0d0 0.0d0 1.d-10 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.07716049382715791d0 2.216394969317974d-12)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGS (INTEGRATION-TEST-F1 2.6d0) 1.0d0
+		      0.0d0 0.0d0 1.d-10 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -5908.755278982137d0 1.2996641844811724d-10)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGS (INTEGRATION-TEST-F11 2.0d0) 1.0d0
+		      1000.0d0 1.d-7 0.0d0 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 5908.755278982137d0 1.2996641844811724d-10)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGS (INTEGRATION-TEST-F11 2.0d0)
+		      1000.0d0 1.0d0 1.d-7 0.0d0 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST -0.3616892186127036d0 3.0167169113304304d-6)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGIU 'INTEGRATION-TEST-F455 0.0d0 0.0d0
+		       0.001d0 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 65536.00000000025d0 7.121667184124547d-4)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGIU (INTEGRATION-TEST-F15 5.0d0) 0.0d0
+		       0.0d0 1.d-7 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0000000000067134d-4 3.0840620210883837d-9)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGIU (INTEGRATION-TEST-F16 1.0d0) 99.9d0
+		       1.d-7 0.0d0 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 2.2758757944687478d0 7.436490135715555d-9)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGI 'INTEGRATION-TEST-MYFN1 1.d-7 0.0d0
+		      1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 2.718281828459045d0 1.58818544083484d-10)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGIL (INTEGRATION-TEST-MYFN2 1.0d0)
+		       1.0d0 1.d-7 0.0d0 1000)))
+  (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+   (LIST 52.740806116727164d0 1.7557038486870624d-4)
+   (MULTIPLE-VALUE-LIST
+    (INTEGRATION-QAGP 'INTEGRATION-TEST-F454
+		      (COPY
+		       (VECTOR 0.0d0 1.0d0 (SQRT 2.0d0)
+			       3.0d0)
+		       'VECTOR-DOUBLE-FLOAT)
+		      0.0d0 0.001d0 1000)))
+  (let ((lisp-unit:*epsilon* 1.0d-9))
+    (LISP-UNIT:ASSERT-NUMERICAL-EQUAL
+     (LIST -0.08994400695837003d0 1.18529017636488d-6)
+     (MULTIPLE-VALUE-LIST
+      (INTEGRATION-QAWC 'INTEGRATION-TEST-F459 -1.0d0 5.0d0
+			0.0d0 0.0d0 0.001d0 1000)))))