Someone asked me this question over email, and I decided to post my answer here since it might be of more general interest.

Suppose that we are trying to compute \int_0^\infty (e^y - 1) \frac{e^{-2.5y}}{y^{1.2}} dy. This integral should be well defined since it decays exponentially as y \to \infty and blows up as 1/y^{0.2} (an integrable singularity) as y \to 0. However, if we naively plug this into QuadGK, it gives a `NaN`

error:

```
julia> quadgk(y -> (exp(y) - 1) * exp(-2.5y) / y^1.2, 0, Inf)
ERROR: DomainError with 0.875:
integrand produced NaN in the interval (0.75, 1.0)
```

Is this a bug in QuadGK? No. In fact, if we evaluate this integrand for large `y`

, it really does give `NaN`

:

```
julia> f(y) = (exp(y)-1) * exp(-2.5y) / y^1.2
f (generic function with 1 method)
julia> f(1000)
NaN
```

Is this a bug in Julia? No. In fact, exactly the same thing will happen in *any* language (Python, Matlab, R, etc.) that uses the usual double-precision floating-point (`Float64`

) arithmetic, due to underflow and overflow. Let’s look at this example more closely. The problem is that `exp(1000) ≈ 1.97e434`

is larger than the largest representable `Float64`

value, so it overflows to `Inf`

, while `exp(-2.5 * 1000) ≈ 1.84e-1086`

is smaller than the smallest representable magnitude so it underflows to `0.0`

:

```
julia> exp(1000) # overflows
Inf
julia> exp(-2.5 * 1000) # underflows
0.0
```

and `Inf * 0.0`

gives `NaN`

— there is no way to evaluate it to a number because the computer has lost track of how large and small the magnitudes are.

How can we avoid this? The usual solution is to refactor your expression to combine the values analytically *before* they over/underflow. In this case, we want to combine e^y e^{-2.5y} = e^{-1.5y}.

However, we should also be careful at y \to 0, since if you compute e^y - 1 naively by `exp(y) - 1`

, it is subject to severe cancellation error for small `y`

— this is why most numerics libraries provide an `expm1`

function to compute it more accurately. So, I would do something like:

```
julia> f(y) = (y < 10 ? expm1(y) * exp(-2.5y) : (exp(-1.5y) - exp(-2.5y))) / y^1.2
f (generic function with 1 method)
julia> quadgk(f, 0, Inf) # default tolerance rtol = √ε ≈ 1e-8
(0.6790524782732376, 9.565804913317641e-9)
julia> quadgk(f, 0, Inf, rtol=1e-14)
(0.6790524809899963, 5.8288819115959285e-15)
```

which uses `expm1`

for sufficiently small arguments and explicitly combines the `exp`

factors for large arguments. The `quadgk`

function now works well even if you require a very small error tolerance close to machine precision.