Quadrature results inconsistent with scipy.integrate.quad and mpmath.quad

[theDDA]

I have a highly oscillatory function that I'd like to integrate on an interval. The function has been omitted as it's quite complicated, but it's a single variable vectorized function.

Running

import quadpy
import mpmath as mp
from scipy.integrate import quad

print("Scipy:")
print(quad(i,0,6,limit=5000,complex_func=True))
print("Mpmath:")
print(mp.fp.quad(i,np.linspace(0,6,1500), error=True))
print("Quadpy:")
print(quadpy.quad(i,0,6, limit=5000))

returns

Scipy:
((-257085873.20190012+274878176.58262545j), (3.788039545497331+4.050770505976548j))
Mpmath:
(array([-2.57085873e+08+2.74878177e+08j]), 0.016388469075099873)
Quadpy:
(np.complex128(162736641.65176928+500993494.40351963j), np.float64(0.6246959267121839))

Scipy and mpmath agree with each other quite well, but I have pretty much no idea how to make quadpy work. I'm not even sure if it's a bug, it may well be that both scipy and mpmath are wrong. All I know is that on every subinterval they agree, while quadpy doesn't seem to be consistent with itself between different schemes.

I've tried quadpy.c1.integrate_adaptive(i,(0,6),minimum_interval_length=6/100000) and combinations of various schemes and degrees. I know that Gauss-Legendre quadrature is used for this function, so I've spent most of my time with that but to no avail.

All of the above holds if I reduce the endpoint of the interval to 6/100 or 6/1000, where the function is much less oscillatory.

I realize this isn't reproducible without the original function, but I was hoping there are some options that I'm missing.

[nschloe]

What version of quadpy are you using?

[theDDA]

It's the conda package, 0.16.10.

FWIW it fails even for intervals where the function appears to be well-behaved.

[nschloe]

The current version is 0.17.21, your issue might be fixed there. We'll have to see about getting conda on track again.