Improper Integrals & Special Functions

If $f$ is Riemann integrable on $[a, b]$ for every $b > a$, we define

$\int_a^{\infty} f(x)\,dx = \lim_{b \to \infty} \int_a^b f(x)\,dx,$

provided the limit exists. If the limit exists and is finite, the improper integral converges; otherwise it diverges.

Similarly, $\int_{-\infty}^b f(x)\,dx = \lim_{a \to -\infty} \int_a^b f(x)\,dx$, and for doubly infinite integrals:

$\int_{-\infty}^{\infty} f(x)\,dx = \int_{-\infty}^c f(x)\,dx + \int_c^{\infty} f(x)\,dx$

for any $c \in \mathbb{R}$, provided both pieces converge. The choice of $c$ doesn’t matter when both parts converge — this follows from the additivity of the Riemann integral over adjacent intervals.

The integral $\int_1^{\infty} \frac{1}{x^p}\,dx$ is the fundamental benchmark for convergence at infinity.

For $p \neq 1$:

$\int_1^b \frac{1}{x^p}\,dx = \frac{x^{1-p}}{1-p}\bigg|_1^b = \frac{b^{1-p} - 1}{1-p}.$

• If $p > 1$: the exponent $1 - p < 0$, so $b^{1-p} \to 0$ as $b \to \infty$. The integral converges to $\frac{1}{p-1}$.
• If $p < 1$: the exponent $1 - p > 0$, so $b^{1-p} \to \infty$. The integral diverges.

For $p = 1$: $\int_1^b \frac{1}{x}\,dx = \ln b \to \infty$. Diverges.

The integral $\int_1^{\infty} \frac{1}{x^p}\,dx$ converges if and only if $p > 1$. This is the $p$-test — every other Type I convergence question reduces to comparing against $1/x^p$.

$\int_0^{\infty} e^{-x}\,dx = \lim_{b \to \infty} \left[-e^{-x}\right]_0^b = \lim_{b \to \infty} (1 - e^{-b}) = 1.$

Exponential decay is more than fast enough to make the area finite. In fact, $e^{-x}$ decays faster than any power $1/x^p$, which is why exponential tails are so well-behaved in probability.

$\int_1^{\infty} \frac{1}{x}\,dx = \lim_{b \to \infty} \ln b = \infty.$

The harmonic function $1/x$ decays, but too slowly — the area accumulates without bound. This is the borderline case $p = 1$ of the $p$-test, and it’s the continuous analog of the divergent harmonic series $\sum 1/n$.

For $\int_{-\infty}^{\infty} f(x)\,dx$, we split at any finite $c$ and require both $\int_{-\infty}^c f$ and $\int_c^{\infty} f$ to converge independently. The Cauchy principal value $\text{PV}\int_{-\infty}^{\infty} f(x)\,dx = \lim_{b \to \infty} \int_{-b}^{b} f(x)\,dx$ is a weaker notion — it can exist even when the improper integral diverges. For example, $\int_{-\infty}^{\infty} x\,dx$ has Cauchy principal value $0$ (by symmetry: $\int_{-b}^{b} x\,dx = 0$ for all $b$), but the improper integral diverges because $\int_0^{\infty} x\,dx = \infty$.

If $f$ is Riemann integrable on $[\varepsilon, b]$ for every $\varepsilon \in (a, b]$ but $f$ is unbounded near $a$, we define

$\int_a^b f(x)\,dx = \lim_{\varepsilon \to a^+} \int_{\varepsilon}^b f(x)\,dx,$

provided the limit exists. Similarly, if $f$ is unbounded near $b$:

$\int_a^b f(x)\,dx = \lim_{\varepsilon \to b^-} \int_a^{b - \varepsilon} f(x)\,dx.$

For singularities at an interior point $c \in (a, b)$, we split: $\int_a^b f = \int_a^c f + \int_c^b f$, and both limits must exist independently.

The integral $\int_0^1 \frac{1}{x^p}\,dx$ is the benchmark for convergence near a singularity.

For $p \neq 1$:

$\int_\varepsilon^1 \frac{1}{x^p}\,dx = \frac{x^{1-p}}{1-p}\bigg|_\varepsilon^1 = \frac{1 - \varepsilon^{1-p}}{1-p}.$

• If $p < 1$: $\varepsilon^{1-p} \to 0$ as $\varepsilon \to 0^+$ (since $1 - p > 0$). The integral converges to $\frac{1}{1-p}$.
• If $p > 1$: $\varepsilon^{1-p} \to \infty$. The integral diverges.

For $p = 1$: $\int_\varepsilon^1 \frac{1}{x}\,dx = -\ln \varepsilon \to \infty$. Diverges.

The integral $\int_0^1 \frac{1}{x^p}\,dx$ converges if and only if $p < 1$. Note the complementary relationship with the Type I $p$-test: $p > 1$ for convergence at infinity, $p < 1$ for convergence at zero. The borderline $p = 1$ diverges in both cases.

$\int_0^1 \frac{1}{\sqrt{x}}\,dx$: here $p = 1/2 < 1$, so the integral converges. Explicitly:

$\int_\varepsilon^1 \frac{1}{\sqrt{x}}\,dx = 2\sqrt{x}\Big|_\varepsilon^1 = 2(1 - \sqrt{\varepsilon}) \to 2 \text{ as } \varepsilon \to 0^+.$

The area under $1/\sqrt{x}$ on $(0, 1]$ is finite ($= 2$) despite the function blowing up at $0$.

If $f$ has a singularity at $c \in (a, b)$, we split: $\int_a^b f = \int_a^c f + \int_c^b f$, and both parts must converge independently. For example:

$\int_{-1}^{1} \frac{1}{|x|^{2/3}}\,dx = \int_{-1}^{0} \frac{1}{|x|^{2/3}}\,dx + \int_0^1 \frac{1}{x^{2/3}}\,dx.$

Each piece converges ($p = 2/3 < 1$), so the integral converges. The total value is $2 \cdot 3 = 6$.

Some integrals involve both an unbounded interval and an unbounded integrand. For example, $\int_0^{\infty} \frac{1}{\sqrt{x}(1+x)}\,dx$ has a Type II singularity at $0$ (where $1/\sqrt{x}$ blows up) and requires a Type I analysis as $x \to \infty$ (where the integrand decays like $1/x^{3/2}$). Split at $x = 1$ and handle each piece separately: near $0$, compare with $1/\sqrt{x}$ ($p = 1/2 < 1$, converges); at infinity, compare with $1/x^{3/2}$ ($p = 3/2 > 1$, converges). The integral converges and equals $\pi$.

Suppose $0 \le f(x) \le g(x)$ for all $x \ge a$.

(a) If $\int_a^{\infty} g(x)\,dx$ converges, then $\int_a^{\infty} f(x)\,dx$ converges, and $\int_a^{\infty} f \le \int_a^{\infty} g$.

(b) If $\int_a^{\infty} f(x)\,dx$ diverges, then $\int_a^{\infty} g(x)\,dx$ diverges.

The analogous result holds for Type II integrals.

For part (a): let $F(b) = \int_a^b f(x)\,dx$ and $G(b) = \int_a^b g(x)\,dx$. By monotonicity of the Riemann integral (Topic 7, Theorem 3), $f \le g$ implies $F(b) \le G(b)$ for all $b > a$. Since $f \ge 0$, the function $F$ is increasing in $b$: adding more non-negative area can only increase the integral. Since $G(b) \to \int_a^\infty g < \infty$, the function $F$ is increasing and bounded above. By the Monotone Convergence principle (Topic 3), $\lim_{b \to \infty} F(b)$ exists and is finite.

Part (b) is the contrapositive of part (a): if $\int g$ converged, then $\int f$ would converge too, contradicting the hypothesis.

$\int_1^{\infty} \frac{1}{1+x^2}\,dx$ converges. For $x \ge 1$:

$\frac{1}{1+x^2} \le \frac{1}{x^2},$

and $\int_1^{\infty} \frac{1}{x^2}\,dx$ converges ($p = 2 > 1$). By the comparison test, $\int_1^{\infty} \frac{1}{1+x^2}\,dx$ converges.

In fact, we can compute the exact value: $\int_1^{\infty} \frac{1}{1+x^2}\,dx = \arctan(\infty) - \arctan(1) = \frac{\pi}{2} - \frac{\pi}{4} = \frac{\pi}{4}$.

Suppose $f(x) > 0$ and $g(x) > 0$ for $x \ge a$, and $\lim_{x \to \infty} \frac{f(x)}{g(x)} = L$.

• If $0 < L < \infty$: $\int_a^{\infty} f$ and $\int_a^{\infty} g$ either both converge or both diverge.
• If $L = 0$ and $\int g$ converges: then $\int f$ converges.
• If $L = \infty$ and $\int g$ diverges: then $\int f$ diverges.

For the case $0 < L < \infty$: choose $\varepsilon = L/2$. There exists $M$ such that for $x \ge M$:

$\frac{L}{2} < \frac{f(x)}{g(x)} < \frac{3L}{2},$

which gives $\frac{L}{2}g(x) < f(x) < \frac{3L}{2}g(x)$ for $x \ge M$. Apply the direct comparison test:

• If $\int g$ converges: $f(x) < \frac{3L}{2}g(x)$, so $\int f$ converges by comparison with $\frac{3L}{2}\int g$.
• If $\int g$ diverges: $f(x) > \frac{L}{2}g(x)$, so $\int f \ge \frac{L}{2}\int g$ diverges.

The cases $L = 0$ and $L = \infty$ follow by similar comparison arguments using one-sided bounds.

$\int_1^{\infty} \frac{x}{x^3 + 1}\,dx$ converges. Compare with $g(x) = \frac{1}{x^2}$:

$\frac{f(x)}{g(x)} = \frac{x/(x^3 + 1)}{1/x^2} = \frac{x^3}{x^3 + 1} \to 1 \text{ as } x \to \infty.$

Since $L = 1 \in (0, \infty)$ and $\int_1^{\infty} \frac{1}{x^2}\,dx$ converges, the limit comparison test gives convergence. We didn’t need to compute the integral — knowing its asymptotic behavior was enough.

An improper integral $\int_a^{\infty} f(x)\,dx$ converges absolutely if $\int_a^{\infty} |f(x)|\,dx$ converges. It converges conditionally if $\int_a^{\infty} f(x)\,dx$ converges but $\int_a^{\infty} |f(x)|\,dx$ diverges.

If $\int_a^{\infty} |f(x)|\,dx$ converges, then $\int_a^{\infty} f(x)\,dx$ converges, and

$\left|\int_a^{\infty} f(x)\,dx\right| \le \int_a^{\infty} |f(x)|\,dx.$

Write $f = f^+ - f^-$ where $f^+(x) = \max(f(x), 0)$ and $f^-(x) = \max(-f(x), 0)$. Then $|f| = f^+ + f^-$, and $0 \le f^+ \le |f|$, $0 \le f^- \le |f|$. Since $\int |f|$ converges, both $\int f^+$ and $\int f^-$ converge by the comparison test (Theorem 1). Then $\int f = \int f^+ - \int f^-$ converges.

The inequality follows: $\left|\int f\right| = \left|\int f^+ - \int f^-\right| \le \int f^+ + \int f^- = \int |f|$.

The Dirichlet integral $\int_0^{\infty} \frac{\sin x}{x}\,dx = \frac{\pi}{2}$ is a famous result (proved via contour integration or Laplace transforms — we state it without proof here). But $\int_0^{\infty} \frac{|\sin x|}{x}\,dx$ diverges: on each interval $[n\pi, (n+1)\pi]$,

$\int_{n\pi}^{(n+1)\pi} \frac{|\sin x|}{x}\,dx \ge \frac{1}{(n+1)\pi}\int_{n\pi}^{(n+1)\pi} |\sin x|\,dx = \frac{2}{(n+1)\pi},$

and $\sum_{n=0}^{\infty} \frac{2}{(n+1)\pi}$ diverges (harmonic series). The integral converges conditionally but not absolutely — positive and negative oscillations cancel just enough to produce a finite result, but the total variation is infinite.

The comparison and limit comparison tests apply to Type II improper integrals with the obvious modifications: compare behavior near the singularity instead of at infinity. For $\int_0^1 f(x)\,dx$ with a singularity at $0$, compare $f(x)$ to $1/x^p$ as $x \to 0^+$; the $p$-test ($p < 1$ for convergence) provides the benchmark.

For $s > 0$, define

$\Gamma(s) = \int_0^{\infty} t^{s-1} e^{-t}\,dt.$

This is a doubly improper integral: the factor $t^{s-1}$ may blow up at $t = 0$ (a Type II singularity when $s < 1$), and the integral extends to $t = \infty$ (Type I).

$\Gamma(s)$ converges for all $s > 0$.

Near $t = 0$: $t^{s-1}e^{-t} \le t^{s-1}$ (since $e^{-t} \le 1$ for $t \ge 0$), and $\int_0^1 t^{s-1}\,dt$ converges if and only if $s - 1 > -1$, i.e., $s > 0$. This is the Type II $p$-test with $p = 1 - s < 1$.

Near $t = \infty$: For any $s > 0$, the exponential decay $e^{-t}$ dominates the polynomial growth $t^{s-1}$. Specifically, for $t$ large enough, $t^{s-1}e^{-t} \le e^{-t/2}$ (since $t^{s-1} \le e^{t/2}$ for $t$ sufficiently large), and $\int_1^{\infty} e^{-t/2}\,dt = 2e^{-1/2}$ converges.

Split the integral at $t = 1$: $\Gamma(s) = \int_0^1 t^{s-1}e^{-t}\,dt + \int_1^{\infty} t^{s-1}e^{-t}\,dt$.

For the first piece: $0 \le t^{s-1}e^{-t} \le t^{s-1}$ on $[0, 1]$ (since $e^{-t} \le 1$). The integral $\int_0^1 t^{s-1}\,dt = \frac{t^s}{s}\big|_0^1 = \frac{1}{s}$ converges for $s > 0$.

For the second piece: we need $t^{s-1}e^{-t/2} \to 0$ as $t \to \infty$. Since $\lim_{t \to \infty} t^{s-1}/e^{t/2} = 0$ for any $s$ (exponential growth dominates polynomial), there exists $T > 1$ such that $t^{s-1} \le e^{t/2}$ for $t \ge T$. Then $t^{s-1}e^{-t} \le e^{-t/2}$ for $t \ge T$, and $\int_T^{\infty} e^{-t/2}\,dt = 2e^{-T/2} < \infty$. The integral $\int_1^T t^{s-1}e^{-t}\,dt$ is a proper Riemann integral (finite interval, bounded integrand), so it converges trivially.

For $s > 0$: $\Gamma(s+1) = s\,\Gamma(s)$.

Integration by parts with $u = t^s$ and $dv = e^{-t}\,dt$:

$\Gamma(s+1) = \int_0^{\infty} t^s e^{-t}\,dt = \left[-t^s e^{-t}\right]_0^{\infty} + s\int_0^{\infty} t^{s-1} e^{-t}\,dt.$

The boundary term vanishes at both limits:

• At $t = 0$: $t^s e^{-t} \to 0$ (for $s > 0$).
• At $t = \infty$: $t^s e^{-t} \to 0$ (exponential decay dominates polynomial growth).

Therefore $\Gamma(s+1) = 0 + s\,\Gamma(s) = s\,\Gamma(s)$.

For positive integers $n$: $\Gamma(n+1) = n!$. Also, $\Gamma(1) = 1$ and $\Gamma(1/2) = \sqrt{\pi}$.

Base case: $\Gamma(1) = \int_0^{\infty} e^{-t}\,dt = 1$.

Induction: By the functional equation, $\Gamma(2) = 1 \cdot \Gamma(1) = 1 = 1!$, $\Gamma(3) = 2 \cdot \Gamma(2) = 2 = 2!$, and generally $\Gamma(n+1) = n \cdot \Gamma(n) = n \cdot (n-1)! = n!$.

Half-integer value: $\Gamma(1/2) = \int_0^{\infty} t^{-1/2} e^{-t}\,dt$. Substitute $t = u^2$, $dt = 2u\,du$:

$\Gamma(1/2) = \int_0^{\infty} (u^2)^{-1/2} e^{-u^2} \cdot 2u\,du = 2\int_0^{\infty} e^{-u^2}\,du = \sqrt{\pi},$

where the last equality is the Gaussian integral (Theorem 7 below).

Using the functional equation $\Gamma(s+1) = s\Gamma(s)$ repeatedly:

$\Gamma(3/2) = \frac{1}{2}\Gamma(1/2) = \frac{\sqrt{\pi}}{2}, \qquad \Gamma(5/2) = \frac{3}{2}\Gamma(3/2) = \frac{3\sqrt{\pi}}{4}.$

In general, $\Gamma(n + 1/2) = \frac{(2n)!}{4^n n!}\sqrt{\pi}$. These half-integer Gamma values appear in the surface area of spheres, the volume of balls in $\mathbb{R}^n$, and the normalizing constants of the Student-$t$ and Chi-squared distributions.

The Gamma function is not just a factorial extension — it’s the factorial extension with the best analytic properties. The Bohr-Mollerup theorem states that $\Gamma$ is the unique function on $(0, \infty)$ satisfying: (i) $\Gamma(1) = 1$, (ii) $\Gamma(s+1) = s\,\Gamma(s)$, and (iii) $\log \Gamma$ is convex. The log-convexity condition rules out other interpolations like $f(s) = \Gamma(s)(1 + \varepsilon\sin(2\pi s))$ that also satisfy (i) and (ii) but oscillate between the integer values.

For $a, b > 0$, define

$B(a, b) = \int_0^1 t^{a-1}(1-t)^{b-1}\,dt.$

This is a doubly improper integral when $a < 1$ or $b < 1$ (Type II singularities at $t = 0$ or $t = 1$). The integral converges for all $a, b > 0$ by the Type II $p$-test at each endpoint: near $t = 0$, the integrand behaves like $t^{a-1}$, and $\int_0^{1/2} t^{a-1}\,dt$ converges iff $a > 0$; near $t = 1$, it behaves like $(1-t)^{b-1}$, and $\int_{1/2}^1 (1-t)^{b-1}\,dt$ converges iff $b > 0$.

For all $a, b > 0$:

$B(a, b) = \frac{\Gamma(a)\,\Gamma(b)}{\Gamma(a+b)}.$

Compute the product $\Gamma(a)\,\Gamma(b)$ as a double integral:

$\Gamma(a)\,\Gamma(b) = \int_0^{\infty} \int_0^{\infty} s^{a-1} t^{b-1} e^{-(s+t)}\,ds\,dt.$

Substitute $s = uv$, $t = u(1-v)$ with $u \in (0, \infty)$ and $v \in (0, 1)$. The Jacobian of this transformation is $\left|\frac{\partial(s,t)}{\partial(u,v)}\right| = u$, so:

$\Gamma(a)\,\Gamma(b) = \int_0^{\infty} \int_0^1 (uv)^{a-1}\bigl(u(1-v)\bigr)^{b-1} e^{-u} \cdot u\,dv\,du.$

Separating:

$= \underbrace{\int_0^{\infty} u^{a+b-1} e^{-u}\,du}_{\Gamma(a+b)} \cdot \underbrace{\int_0^1 v^{a-1}(1-v)^{b-1}\,dv}_{B(a,b)}.$

Therefore $\Gamma(a)\,\Gamma(b) = \Gamma(a+b) \cdot B(a,b)$, which gives $B(a,b) = \frac{\Gamma(a)\,\Gamma(b)}{\Gamma(a+b)}$.

(This proof uses Fubini’s theorem and a 2D change of variables — multivariable techniques that will be developed rigorously in Track 4. We preview them here because the result is too important to postpone.)

$B(1/2, 1/2) = \frac{\Gamma(1/2)^2}{\Gamma(1)} = \frac{(\sqrt{\pi})^2}{1} = \pi.$

Direct verification: $\int_0^1 \frac{1}{\sqrt{t(1-t)}}\,dt = \int_0^1 \frac{dt}{\sqrt{t - t^2}}$. Completing the square and substituting $t = \frac{1}{2}(1 + \sin\theta)$ gives $\pi$.

$B(a, 1) = \int_0^1 t^{a-1}\,dt = \frac{1}{a}$.

Verify via the Gamma relationship: $\frac{\Gamma(a)\,\Gamma(1)}{\Gamma(a+1)} = \frac{\Gamma(a) \cdot 1}{a\,\Gamma(a)} = \frac{1}{a}$. Similarly, $B(1, b) = 1/b$.

If $X \sim \text{Beta}(\alpha, \beta)$, its density is $f(x) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)}$ for $x \in (0, 1)$. The Beta function is precisely the normalizing constant that makes $\int_0^1 f(x)\,dx = 1$. The parameters $\alpha, \beta > 0$ control the shape: $\alpha = \beta = 1$ gives the uniform distribution on $[0, 1]$; $\alpha = \beta \gg 1$ concentrates the density around $x = 1/2$; $\alpha > \beta$ skews the density toward $x = 1$.

In Bayesian inference, the Beta distribution is the conjugate prior for the Binomial likelihood. If you observe $k$ successes in $n$ trials and your prior is $\text{Beta}(\alpha, \beta)$, the posterior is $\text{Beta}(\alpha + k, \beta + n - k)$. The marginal likelihood involves the ratio $B(\alpha + k, \beta + n - k) / B(\alpha, \beta)$, which reduces to a ratio of Gamma functions.

$\int_{-\infty}^{\infty} e^{-x^2}\,dx = \sqrt{\pi}.$

Equivalently, $\int_0^{\infty} e^{-x^2}\,dx = \frac{\sqrt{\pi}}{2}$.

Let $I = \int_0^{\infty} e^{-x^2}\,dx$. The key trick is to compute $I^2$ as a double integral:

$I^2 = \left(\int_0^{\infty} e^{-x^2}\,dx\right)\left(\int_0^{\infty} e^{-y^2}\,dy\right) = \int_0^{\infty}\int_0^{\infty} e^{-(x^2 + y^2)}\,dx\,dy.$

Convert to polar coordinates: $x = r\cos\theta$, $y = r\sin\theta$, with Jacobian $r$, and $x^2 + y^2 = r^2$:

$I^2 = \int_0^{\pi/2}\int_0^{\infty} e^{-r^2}\,r\,dr\,d\theta = \frac{\pi}{2}\int_0^{\infty} r\,e^{-r^2}\,dr.$

The inner integral evaluates by substituting $u = r^2$, $du = 2r\,dr$:

$\int_0^{\infty} r\,e^{-r^2}\,dr = \frac{1}{2}\int_0^{\infty} e^{-u}\,du = \frac{1}{2}.$

Therefore $I^2 = \frac{\pi}{2} \cdot \frac{1}{2} = \frac{\pi}{4}$, so $I = \frac{\sqrt{\pi}}{2}$, and $\int_{-\infty}^{\infty} e^{-x^2}\,dx = 2I = \sqrt{\pi}$.

(This proof uses a double integral and polar coordinates — multivariable tools that will be developed rigorously in Track 4. We preview them here because the result is too important and the proof too elegant to postpone. The reader can follow the geometry: squaring the integral converts a 1D problem into a 2D problem where the radial symmetry of $e^{-(x^2+y^2)}$ makes the integral tractable.)

(a) For $a > 0$: $\int_{-\infty}^{\infty} e^{-ax^2}\,dx = \sqrt{\pi/a}$. (Substitute $u = \sqrt{a}\,x$.)

(b) Completing the square: $\int_{-\infty}^{\infty} e^{-ax^2 + bx}\,dx = \sqrt{\pi/a}\,e^{b^2/(4a)}$.

(c) Even moments: $\int_{-\infty}^{\infty} x^{2n} e^{-x^2}\,dx = \frac{(2n)!}{4^n n!}\sqrt{\pi}$.

These follow from (a) by differentiation with respect to $a$ (for the moments) or completing the square in the exponent.

The Gaussian density is $f(x) = \frac{1}{\sigma\sqrt{2\pi}}\,e^{-(x-\mu)^2/(2\sigma^2)}$. We verify $\int_{-\infty}^{\infty} f(x)\,dx = 1$:

Substitute $u = (x - \mu)/(\sigma\sqrt{2})$, so $dx = \sigma\sqrt{2}\,du$:

$\int_{-\infty}^{\infty} f(x)\,dx = \frac{1}{\sigma\sqrt{2\pi}} \cdot \sigma\sqrt{2} \int_{-\infty}^{\infty} e^{-u^2}\,du = \frac{\sqrt{2}}{\sqrt{2\pi}} \cdot \sqrt{\pi} = 1. \quad \checkmark$

The $\sqrt{2\pi}$ in the denominator of the Gaussian density is there precisely because $\int e^{-x^2}\,dx = \sqrt{\pi}$ — it’s the Gaussian integral doing the normalizing.

The error function $\text{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt$ is the normalized incomplete Gaussian integral. The Gaussian CDF is:

$\Phi(x) = \frac{1}{2}\left[1 + \text{erf}\left(\frac{x}{\sqrt{2}}\right)\right].$

The error function has no closed-form antiderivative — $e^{-t^2}$ cannot be integrated in terms of elementary functions. But it can be computed numerically to machine precision by standard libraries (scipy.special.erf, torch.special.erf). The rapid convergence of the Gaussian tail ($\Phi(x) \approx 1 - \frac{1}{x\sqrt{2\pi}}e^{-x^2/2}$ for large $x$, known as Mill’s ratio) is what makes sub-Gaussian concentration inequalities so powerful.

$n! \sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n \quad \text{as } n \to \infty.$

More precisely: $\lim_{n \to \infty} \frac{n!}{\sqrt{2\pi n}(n/e)^n} = 1$. The relative error is $O(1/n)$:

$n! = \sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1 + \frac{1}{12n} + O(1/n^2)\right).$

We sketch the proof via the Gamma function and Laplace’s method. Write $n! = \Gamma(n+1) = \int_0^{\infty} t^n e^{-t}\,dt$. The integrand $h(t) = t^n e^{-t}$ has its maximum at $t = n$ (set $h'(t) = 0$: $nt^{n-1}e^{-t} - t^n e^{-t} = 0$ gives $t = n$).

Substitute $t = n + \sqrt{n}\,u$ to center the integrand at its maximum. Taylor-expanding $\log h(t)$ around $t = n$:

$\log h(n + \sqrt{n}\,u) = \log h(n) + 0 \cdot (\sqrt{n}\,u) - \frac{u^2}{2} + O(u^3/\sqrt{n}).$

So $h(n + \sqrt{n}\,u) \approx h(n) \cdot e^{-u^2/2}$ for the leading term. Integrating:

$\Gamma(n+1) \approx h(n) \cdot \sqrt{n} \int_{-\infty}^{\infty} e^{-u^2/2}\,du = n^n e^{-n} \cdot \sqrt{n} \cdot \sqrt{2\pi} = \sqrt{2\pi n}\left(\frac{n}{e}\right)^n.$

The Laplace approximation replaces the integrand by a Gaussian centered at its maximum — the $\sqrt{2\pi}$ factor is precisely the Gaussian integral. The correction terms (the $1/(12n)$ and higher) come from including the cubic and higher Taylor terms of $\log h$.

The approximation improves as $n$ increases — the relative error decreases like $1/(12n)$.

In practice, the most useful form is:

$\log(n!) \approx n\log n - n + \frac{1}{2}\log(2\pi n).$

This is the form used throughout information theory. The key application: $\log\binom{n}{k} = \log(n!) - \log(k!) - \log((n-k)!)$. Applying Stirling to each factorial:

$\log\binom{n}{k} \approx nH(k/n),$

where $H(p) = -p\log p - (1-p)\log(1-p)$ is the binary entropy function. Stirling’s approximation converts combinatorial counting problems into entropy calculations — this is the foundation of the “method of types” in information theory.