\Version{GAUSS}                    \Routid{D103}
\Keywords{ADAPTIVE NUMERICAL INTEGRATION GAUSS QUADRATURE}
\Author{<NAME>}             \Library{MATHLIB}
\Submitter{}                       \Submitted{02.05.1966}
\Language{Fortran}                  \Revised{15.03.1993}
\Cernhead{Adaptive Gaussian Quadrature}
Function subprograms {\tt GAUSS}, {\tt DGAUSS} and {\tt QGAUSS} compute,
to an attempted specified accuracy, the value of the integral
 $$ I=\int^B_A f(x)dx.$$
\par
The quadruple-precision version {\tt QGAUSS} is available only on
computers which support a {\tt REAL*16} Fortran data type.
\Structure
{\tt FUNCTION} subprograms  \\
User Entry Names: \Rdef{GAUSS}, \Rdef{DGAUSS}, \Rdef{QGAUSS}\\
Files Referenced: {\tt Unit 6} \\
External References: \Rind{MTLMTR}{N002},  \Rind{ABEND}{Z035},
user-supplied {\tt FUNCTION} subprogram
\Usage
In any arithmetic expression,
\begin{center}
{\tt GAUSS(F,A,B,EPS)}, \quad {\tt DGAUSS(F,A,B,EPS)} \quad or \quad
{\tt QGAUSS(F,A,B,EPS)}
\end{center}
has the approximate value of the integral $I$.
\begin{DLtt}{123456}
\item[F] Name of a user-supplied {\tt FUNCTION} subprogram, declared
{\tt EXTERNAL} in the calling program. This subprogram must set
$\mathtt{F(X)} = f(\mathtt{X})$.
\item [A,B]End-points of integration interval. Note that {\tt B}
may be less than {\tt A}.
\item[EPS]Accuracy parameter (see {\bf Accuracy}).
\end{DLtt}
{\tt GAUSS} is of type {\tt REAL}, {\tt DGAUSS} is of type
{\tt DOUBLE PRECISION}, {\tt QGAUSS} is of type {\tt REAL*16},
and the arguments {\tt F}, {\tt A}, {\tt B},
{\tt EPS} and {\tt X} (in {\tt F}) have the same type as the function
name.
\Method
For any interval $[a,b]$ we define $g_8(a,b)$ and $g_{16}(a,b)$ to be the
8-point and 16-point Gaussian quadrature approximations to
     $$\int^b_a f(x)dx $$
and define
$$ r(a,b) =\frac{|g_{16} (a,b) - g_8(a,b)|}{1+|g_{16}(a,b)|}. $$
Then, with $G$ = {\tt GAUSS} or {\tt DGAUSS},
$$ G =\sum_{i=1}^kg_{16}(x_{i-1},x_i),$$
\newpage
where, starting with $ x_0=A $ and finishing with $x_k=B$,
the subdivision points $ x_i \ (i=1,2,\ldots) $ are given by
$$ x_i = x_{i-1} + \lambda (B-x_{i-1}), $$
with $\lambda$ equal to the first member of the sequence
$1,\frac{1}{2},\frac{1}{4},\ldots $ for which
$r(x_{i-1},x_i) < \mathtt{EPS} $.
If, at any stage in the process of subdivision, the ratio
  $$ q=\left|\frac{x_i-x_{i-1}}{B-A}\right| $$
is so small that $1+0.005q$ is indistinguishable from 1 to
machine accuracy, an error exit occurs with the function value
set equal to zero.
\Accuracy
Unless there is severe cancellation of positive and negative
values of $f(x)$ over the interval $[A,B]$, the argument {\tt EPS}
may be considered as specifying a bound on the {\it relative} error of
$I$ in the case $ |I|>1$, and a bound on the {\it absolute} error in
the case $|I|<1$. More precisely, if $k$ is the number of sub-intervals
contributing to the approximation (see Method), and if
$$ I_{abs} = \int^B_A|f(x)|dx,$$
then the relation
$$ \frac{|G - I|}{I_{abs}+k}< \mathtt{EPS} $$
will nearly always be true, provided the routine terminates
without printing an error message. For functions
$f$ having no singularities in the closed interval $[A,B]$
the accuracy will usually be much higher than this.
\Errorh
Error {\tt D103.1:}  The requested accuracy cannot be
obtained (see {\bf Method}).
The function value is set equal to zero, and a message is written on
{\tt Unit 6} unless subroutine {\tt MTLSET} (N002) has been called.
\Notes
Values of the function $f(x)$ at the interval end-points
$A$ and $B$ are not required. The subprogram may therefore
be used when these values are undefined.
\\ $\bullet$
