calculus and analysis - Integrate returns imaginary answer for smooth, real function

Bug introduced in 7 or earlier and persisting through 12.0.0.0 or later

---

I'm trying to evaluate the integral:

$$\int_0^{\infty} \frac{1}{4 b \sqrt{\pi} r} e^{-(b-r)^2}(e^{4 b r} - 1) \mathrm{d}r$$

with $b>0$.

Integrate[(E^-(b + r)^2 (-1 + 
      E^(4 b r)))/(4 b Sqrt[π] r), {r, 0, ∞}, 
 Assumptions -> b > 0]
(* -((I E^-b^2 Sqrt[π])/(4 b)) *)

Unfortunately, the function being integrated is purely real and smooth, so the answer should be real, but the one provided by Mathematica is purely imaginary.

-((I E^-b^2 Sqrt[π])/(4 b)) /. b -> 1/(π + E) // N
(* 0. - 2.52206 I *)

NIntegrate is able to correctly return real values.

NIntegrate[-(E^-(b + r)^2/(4 b Sqrt[π] r)) + 
   E^(4 b r - (b + r)^2)/(4 b Sqrt[π] r) /. 
  b -> 1/(π + E), {r, 0, ∞}]
(* 0.490405 *)

Unlike some other cases I have seen, Mathematica does not know the elementary antiderivative for this function, so I cannot tell whether some issue with discontinuities or branch cuts is causing this problem. What explains this strange result; is it a bug?

I've checked that versions 9.0, 10.2, and 10.3 all give the same result.