Compute the integral
b
/
[
I f(x) w(x) dx
]
/
a

where w(x) is the weight function below.

f - function subprogram defining the integrand
function f(x). the actual name for f needs to be
declared e x t e r n a l in the driver program.

a - lower limit of integration

b - upper limit of integration, b.gt.a
if b.le.a, the routine will end with ier = 6.

alfa - parameter in the integrand function, alfa.gt.(-1)
if alfa.le.(-1), the routine will end with
ier = 6.

beta - parameter in the integrand function, beta.gt.(-1)
if beta.le.(-1), the routine will end with
ier = 6.

integr - indicates which weight function is to be used
= 1 (x-a)**alfa*(b-x)**beta
= 2 (x-a)**alfa*(b-x)**beta*log(x-a)
= 3 (x-a)**alfa*(b-x)**beta*log(b-x)
= 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
if integr.lt.1 or integr.gt.4, the routine
will end with ier = 6.

epsabs - absolute accuracy requested
epsrel - relative accuracy requested
if epsabs.le.0
and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
the routine will end with ier = 6.

limit - dimensioning parameter for iwork
limit determines the maximum number of
subintervals in the partition of the given
integration interval (a,b), limit.ge.2.
if limit.lt.2, the routine will end with ier = 6.
RETURN
result - double precision
approximation to the integral

ier - integer
ier = 0 normal and reliable termination of the
routine. it is assumed that the requested
accuracy has been achieved.
ier.gt.0 abnormal termination of the routine
the estimates for the integral and error
are less reliable. it is assumed that the
requested accuracy has not been achieved.
abserr - double precision
estimate of the modulus of the absolute error,
which should equal or exceed abs(i-result)

neval - integer
number of integrand evaluations

last - on return, last equals the number of
subintervals produced in the subdivision process,
which determines the significant number of
elements actually in the work arrays.

ERROR MESSAGES
ier = 1 maximum number of subdivisions allowed
has been achieved. one can allow more
subdivisions by increasing the value of
limit (and taking the according dimension
adjustments into account). however, if
this yields no improvement it is advised
to analyze the integrand, in order to
determine the integration difficulties
which prevent the requested tolerance from
being achieved. in case of a jump
discontinuity or a local singularity
of algebraico-logarithmic type at one or
more interior points of the integration
range, one should proceed by splitting up
the interval at these points and calling
the integrator on the subranges.
= 2 the occurrence of roundoff error is
detected, which prevents the requested
tolerance from being achieved.
= 3 extremely bad integrand behaviour occurs
at some points of the integration
interval.
= 6 the input is invalid, because
b.le.a or alfa.le.(-1) or beta.le.(-1) or
or integr.lt.1 or integr.gt.4 or
(epsabs.le.0 and
epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
or limit.lt.2 or lenw.lt.limit*4.
result, abserr, neval, last are set to
zero. except when lenw or limit is invalid
iwork(1), work(limit*2+1) and
work(limit*3+1) are set to zero, work(1)
is set to a and work(limit+1) to b.