Cauchy principal value integrals

[fmeirinhos]

Here's an implementation that can evaluate Cauchy principal value integrals using the method presented here. It does not support infinite boundaries due to this already known problem.

I tried to stick to / reuse most of the base code of do_quadgk but some fundamental changes were needed to account for the evaluation of a function with multiple simple poles. One could argue that that is a bit redundant, since the user could just split the the integral into several 1-simple-poled functions, but I decided to keep it anyway.

I am sure certain things could be polished and a downside is that uses the FastTransforms library for this computation.

Some parts could still be polished (such as supporting functions beyond C -> C) but I think I should hear your feedback regarding bigger changes one should worry about first

[giordano]

a downside is that uses the FastTransforms

I'm not sure this is desirable. This package is purely written in Julia and has minimal dependencies -- only one non-stdlib package so far. FastTransforms.jl instead is larger and has dependencies on multiple binary libraries -- like OpenSpecfun_jll, FFTW_jll, FastTransforms_jll.

[fmeirinhos]

a downside is that uses the FastTransforms

I'm not sure this is desirable.

I have now removed the need for FastTransforms of the Cauchy integrator.

The code seems to be significantly slower (~50%) since the Chebyshev coefficients are calculated with a O(N^2) formula instead of the O(N log N) discrete cosine transform algorithm used by FastTransforms (calling FFTW)

[stevengj]

Why can't one use a transformation like: in which the integral on the right-hand-side has no singularity (assuming f is Lipshitz continuous at x=0) and hence can be done normally via quadgk. (That is, can't you just add and subtract f(0)/x from the integrand?) Am I forgetting something?

(You have to be a little careful about evaluating the integrand at x=0, since you don't want to compute 0/0. Either take the limit, or simply arrange the quadrature rule so that you don't evaluate exactly that point. For example, in the common case where a=-b, you could integrate (f(x) - f(-x))/x from 0 to b, since QuadGK never evaluates the integrand exactly at the endpoints.)

[fmeirinhos]

That is an excellent transformation which I unfortunately did not know about.

I still had the hope that the modified Clenshaw-Curtis quadrature method could outperform quadgk for nastier single-poled integrals (such that the work wouldn't go to waste) but quadgk is rightfully resilient and fast.

May I suggest adding the transformation to the documentation, possibly with a snippet such as

function quadgk_cauchy(f, a, c, b)
  fc = f(c)
  g(x) = (f(x) - fc) / (x - c)
  return quadgk(g, a, c, b)[1] + fc * log(abs((b - c)/(a - c)))
end