title: math range error when trying to integrate a function using scipy quad and solving using brentq
body: <p>I'd like scipy to calculate following integral for me:</p>
  <pre><code>def F(x, theta, mu, sigma, r):
      #equation 1
      def integrand(u):
          return u**(r/mu - 1) * exp(sqrt(2*mu / sigma**2) * (x-theta)*u - u**2/2)
      return scipy.integrate.quad(integrand, 1e-5, np.inf)[0]
  
  def Prime(f, x, theta, mu, sigma, r, h=1e-5):
      # given f, estimates f'(x) using the difference quotient formula 
      return (f(x+h, theta, mu, sigma, r) - f(x, theta, mu, sigma, r)) / h
  
  def b_star(theta, mu, sigma, r, c):
      #equation 2
      def func(b):
          return F(b, theta, mu, sigma, r) - (b-c)*Prime(F, b, theta, mu, sigma, r)
      
      return scipy.optimize.brentq(func, -2, 3)
  </code></pre>
  <p>which correspond with the following equations
  <img src="https://i.stack.imgur.com/XNJfb.png" alt="equation 1" />,</p>
  <p><img src="https://i.stack.imgur.com/gS0LE.png" alt="equation 2" /></p>
  <p>in which I solve equation 2 for b.
  This seems to work most of the time however in some circumstances the following error occurs:</p>
  <pre><code>In: b_star(0.573085,3.213728,0.137655,0,0.001)
  Out:
  ---------------------------------------------------------------------------
  OverflowError                             Traceback (most recent call last)
  &lt;ipython-input-164-1aeb31cb3219&gt; in &lt;module&gt;
  ----&gt; 1 b_star(0.573085,3.213728,0.137655,0,0.001)
  
  &lt;ipython-input-163-14ce7e008b20&gt; in b_star(theta, mu, sigma, r, c)
       38         return F(b, theta, mu, sigma, r) - (b-c)*Prime(F, b, theta, mu, sigma, r)
       39 
  ---&gt; 40     return so.brentq(func, -2, 3)
       41 
       42 def V(x, theta, mu, sigma, r, c):
  
  /Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/scipy/optimize/zeros.py in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
      778     if rtol &lt; _rtol:
      779         raise ValueError(&quot;rtol too small (%g &lt; %g)&quot; % (rtol, _rtol))
  --&gt; 780     r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
      781     return results_c(full_output, r)
      782 
  
  &lt;ipython-input-163-14ce7e008b20&gt; in func(b)
       36 
       37     def func(b):
  ---&gt; 38         return F(b, theta, mu, sigma, r) - (b-c)*Prime(F, b, theta, mu, sigma, r)
       39 
       40     return so.brentq(func, -2, 3)
  
  &lt;ipython-input-163-14ce7e008b20&gt; in F(x, theta, mu, sigma, r)
       18     def integrand(u):
       19         return u**(r/mu - 1) * exp(sqrt(2*mu / sigma**2) * (x-theta)*u - u**2/2)
  ---&gt; 20     return si.quad(integrand, 1e-5, np.inf)[0]
       21 
       22 def G(x, theta, mu, sigma, r):
  
  /Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
      339 
      340     if weight is None:
  --&gt; 341         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
      342                        points)
      343     else:
  
  /Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
      453             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
      454         else:
  --&gt; 455             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
      456     else:
      457         if infbounds != 0:
  
  &lt;ipython-input-163-14ce7e008b20&gt; in integrand(u)
       17     # equation 3.3
       18     def integrand(u):
  ---&gt; 19         return u**(r/mu - 1) * exp(sqrt(2*mu / sigma**2) * (x-theta)*u - u**2/2)
       20     return si.quad(integrand, 1e-5, np.inf)[0]
       21 
  
  OverflowError: math range error
  
  </code></pre>
  <p>I tried several solutions like using the Decimal package or the gmpy package to allow for bigger floating point numbers but the output becomes 0.0 in those cases.
  I also tried converting the inputs to a np.float128 datatype which didn't solve the problem either.</p>
  <p>Is there a fix for this?</p>