Compute the Cauchy principal value of the integral:

b
/
[ f(x)
I ----- dx
] x - c
/
a

f - function subprogram defining the integrand
function f(x).

a - under limit of integration

b - upper limit of integration

c - parameter in the weight function, c.ne.a, c.ne.b.
if c = a or c = b, 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 - limit determines the maximum number of subintervals
in the partition of the given integration interval
(a,b), limit.ge.1.
if limit.lt.1, 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 integral and error are
less reliable. it is assumed that the
requested accuracy has not been achieved.
abserr - double precision
estimate or the modulus of the absolute error,
which should equal or exceed abs(i-result)

neval - integer
number of integrand evaluations

ERROR MESSAGES
ier = 1 maximum number of subdivisions allowed
has been achieved. one can allow more sub-
divisions 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.
if the position of a local difficulty
can be determined (e.g. singularity,
discontinuity within the interval) one
will probably gain from splitting up the
interval at this point and calling
appropriate integrators on the subranges.
= 2 the occurrence of roundoff error is detec-
ted, 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
c = a or c = b or
(epsabs.le.0 and
epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
or limit.lt.1 or lenw.lt.limit*4.
result, abserr, neval, last are set to
zero. exept 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.