I am integrating a function to implement the functions in this paper and I observe different outputs from scipy.integrate.quad and mpmath.quad as per the following:

import mpmath as mp
from scipy import integrate

print(mp.quad(lambda x: x**58 * (1-x)**167 * (x+0.07053941908706418)**86, [0,1]))
print(integrate.quad(lambda x: x**58 * (1-x)**167 * (x+0.07053941908706418)**86, 0, 1))

# Results
# 3.99873917020827e-90
# (1.0583370200366968e-89, 1.5341117982755304e-89)

Based on external considerations, I have reason to believe that the SciPy result is correct (I performed checks on the PDF returned and the values returned by SciPy result is a valid PDF that integrates to 1 over the support).

Am I calling either function, or specifically mp.quad incorrectly?