Understanding Mathematica's formula for $ \int_0^{\infty } x^a \exp \left(-\frac{c x^2+f x}{b}\right) \, dx $ My goal is to integrate the following function:
$$
\int_0^{\infty } x^a \exp \left(-\frac{c x^2+f x}{b}\right) \, dx
$$
where, $a, b, c > 0$ and $a, b, c, f \in \mathbb{R}$.
Mathematica gives me the answer:
$$
\frac{\left(\frac{c}{b}\right)^{-\frac{a}{2}} \left(-f \Gamma \left(1+\frac{a}{2}\right) \, _1F_1\left(1+\frac{a}{2};\frac{3}{2};\frac{f^2}{4 b c}\right)+b \sqrt{\frac{c}{b}} \Gamma \left(\frac{1+a}{2}\right) \,
   _1F_1\left(\frac{1+a}{2};\frac{1}{2};\frac{f^2}{4 b c}\right)\right)}{2 c}
$$
but I don't understand how it obtains this result, or where Kummer's confluent hypergeometric $_1F_1$ comes from. I know some of the integral representations of $_1F_1$, but those range usually integrate from 0 to 1. Can anybody explain how Mathematica solves this? Thanks in advance!

Let $I(a)$ denote your integral. Integration by parts gives
$$
\begin{align*}
I(a) 
& = -\int_{0}^{\infty} \frac{x^{a+1}}{a+1}(-\frac{2cx + f}{b})\exp(-\frac{cx^2 + fx}{b}) \\
& = \frac{2c}{b(a+1)}I(a+2) + \frac{f}{b(a+1)}I(a+1)
\end{align*}
$$
or, in other terms,
$$
I(a+2) = - \frac{f}{2c} I(a+1) + \frac{b(a+1)}{2c}I(a)
$$
I assume that Mathematica tries to solve this recurrence (computing $I(0)$ and $I(1)$ along the way) and ends up with hypergeometric functions because they satisfy many.
