A function with vectorized intervals does not provide correct answers for 1D line segment

[bigm128]

I was trying to compute a function with multiple intervals but it keeps failing or gives wrong answers. Here is an example:

val, err = quadpy.c1.integrate_adaptive(
    lambda x: x,
    [[0.0, 0.0, 0.0], [1.0, 2.0, 3.0]]
)

which returns:

Traceback (most recent call last):
  File "C:\...\test_quadpy.py", line 11, in <module>
    val, err = quadpy.quad(lambda x: x, [0.0, 0.0, 0.0], [1.0, 2.0, 3.0])
  File "C:\...\quadpy\_scipy_compat.py", line 55, in quad
    val, err = integrate_adaptive(
  File "C:\...\quadpy\c1\_adaptive.py", line 141, in integrate_adaptive
    good_values_sum = good_values_sum.reshape(orig_shape)
ValueError: cannot reshape array of size 1 into shape (3,)

But if I specify the domain_shape and range_shape in the function call like this:

val, err = quadpy.c1.integrate_adaptive(
    lambda x: x,
    [[0.0, 0.0, 0.0], [1.0, 2.0, 3.0]],
    domain_shape=(3,),
    range_shape=(3,)
)

the function returns [1.87082869 3.74165739 5.61248608] with error estimates [7.75540290e-20 1.55108058e-19 2.32662087e-19] which are not the correct answers (should be [0.5, 0.333, 1.5]).

I'm using numpy version 1.22.3 and scipy version 1.8.0.

Thank you so much.