Evaluate definite interval of normal distribution I know that an easy to handle formula for the CDF of a normal distribution is somewhat missing, due to the complicated error function in it.
However, I wonder if there is a a nice formula for $N(c_{-} \leq x < c_{+}| \mu, \sigma^2)$. Or what the "state of the art" approximation for this problem might be.
 A: It depends on exactly what you are looking for. Below are some brief details and references.
Much of the literature for approximations centers around the function
$$
Q(x) = \int_x^\infty \frac{1}{\sqrt{2\pi}} e^{-\frac{u^2}{2}} \, \mathrm{d}u
$$
for $x > 0$. This is because the function you provided can be decomposed as a simple difference of the function above (possibly adjusted by a constant). This function is referred to by many names, including "upper-tail of the normal distribution", "right normal integral", and "Gaussian $Q$-function", to name a few. You'll also see approximations to Mills' ratio, which is
$$
R(x) = \frac{Q(x)}{\varphi(x)}
$$
where $\varphi(x) = (2\pi)^{-1/2} e^{-x^2 / 2}$ is the Gaussian pdf.
Here I list some references for various purposes that you might be interested in.
Computational
The de-facto standard for computing the $Q$-function or the related complementary error function is

W. J. Cody, Rational Chebyshev Approximations for the Error Function, Math. Comp., 1969, pp. 631--637.

Every (self-respecting) implementation uses this paper. (MATLAB, R, etc.)
"Simple" Approximations
Abramowitz and Stegun have one based on a polynomial expansion of a transformation of the input. Some people use it as a "high-precision" approximation. I don't like it for that purpose since it behaves badly around zero. For example, their approximation does not yield $\hat{Q}(0) = 1/2$, which I think is a big no-no. Sometimes bad things happen because of this.
Borjesson and Sundberg give a simple approximation which works pretty well for most applications where one only requires a few digits of precision. The absolute relative error is never worse than 1%, which is quite good considering its simplicity. The basic approximation is
$$
\hat{Q}(x) = \frac{1}{(1-a) x + a \sqrt{x^2 + b}} \varphi(x)
$$
and their preferred choices of the constants are $a = 0.339$ and $b = 5.51$. That reference is

P. O. Borjesson and C. E. Sundberg. Simple approximations of the error function Q(x) for communications
applications. IEEE Trans. Commun., COM-27(3):639–643, March 1979.

Here is a plot of its absolute relative error.

The electrical-engineering literature is awash with various such approximations and seem to take an overly intense interest in them. Many of them are poor though or expand to very strange and convoluted expressions.
You might also look at

W. Bryc. A uniform approximation to the right normal integral. Applied Mathematics and Computation,
127(2-3):365–374, April 2002.

Laplace's continued fraction
Laplace has a beautiful continued fraction which yields successive upper and lower bounds for every value of $x > 0$. It is, in terms of Mills' ratio,
$$
R(x) = \frac{1}{x+}\frac{1}{x+}\frac{2}{x+}\frac{3}{x+}\cdots ,
$$
where the notation I've used is fairly standard for a continued fraction, i.e., $1/(x+1/(x+2/(x+3/(x+\cdots))))$. This expression doesn't converge very fast for small $x$, though, and it diverges at $x = 0$.
This continued fraction actually yields many of the "simple" bounds on $Q(x)$ that were "rediscovered" in the mid-to-late 1900s. It's easy to see that for a continued fraction in "standard" form (i.e., composed of positive integer coefficients), truncating the fraction at odd (even) terms gives an upper (lower) bound.
Hence, Laplace tells us immediately that
$$
\frac{x}{x^2 + 1} < R(x) < \frac{1}{x} \>,
$$
both of which are bounds that were "rediscovered" in the mid-1900's. In terms of the $Q$-function, this is equivalent to
$$
\frac{x}{x^2 + 1} \varphi(x) < Q(x) < \frac{1}{x} \varphi(x) .
$$
An alternative proof of this using simple integration by parts can be found in S. Resnick, Adventures in Stochastic Processes, Birkhauser, 1992, in Chapter 6 (Brownian motion). The absolute relative error of these bounds is no worse than $x^{-2}$, as shown in this related answer.
Notice, in particular, that the inequalities above immediately imply that $Q(x) \sim \varphi(x)/x$. This fact can be established using L'Hopital's rule as well. This also helps explain the choice of the functional form of the Borjesson-Sundberg approximation. Any choice of $a \in [0,1]$ maintains the asymptotic equivalence as $x \to \infty$. The parameter $b$ serves as a "continuity correction" near zero.
Here is a plot of the $Q$-function and the two Laplace bounds.

C-I. C. Lee has a paper from the early 1990's that does a "correction" for small values of $x$. See

C-I. C. Lee. On Laplace continued fraction for the normal integral. Ann. Inst. Statist. Math., 44(1):107–120,
March 1992.


Durrett's Probability: Theory and Examples provides the classical upper and lower bounds on $Q(x)$ on pages 6–7 of the 3rd edition. They're meant for larger values of $x$ (say, $x > 3$) and are asymptotically tight.
Hopefully this will get you started. If you have a more specific interest, I might be able to point you somewhere.
A: I suppose I'm too late the hero, but I wanted to comment on cardinal's post, and this comment became too big for its intended box.
For this answer, I'm assuming $x >0$; appropriate reflection formulae can be used for negative $x$.
I'm more used to dealing with the error function $\mathrm{erf}(x)$ myself, but I'll try to recast what I know in terms of Mills's ratio $R(x)$ (as defined in cardinal's answer).
There are in fact alternative ways for computing the (complementary) error function apart from using Chebyshev approximations. Since the use of a Chebyshev approximation requires the storage of not a few coefficients, these methods might have an edge if array structures are a bit costly in your computing environment (you could inline the coefficients, but the resulting code would probably look like a baroque mess).
For "small" $|x|$, Abramowitz and Stegun give a nicely behaved series (at least better behaved than the usual Maclaurin series):
$$R(x)=\sqrt{\frac{\pi}{2}}\exp\left(\frac{x^2}{2}\right)-x\sum_{j=0}^\infty\frac{2^j j!}{(2j+1)!}x^{2j}$$
(adapted from formula 7.1.6)
Note that the coefficients of $x^{2j}$ in the series $c_j=\frac{2^j j!}{(2j+1)!}$ can be computed by starting with $c_0=1$ and then using the recursion formula $c_{j+1}=\frac{c_j}{2j+3}$. This is convenient when implementing the series as a summation loop.

cardinal gave the Laplacian continued fraction as a way to bound Mills's ratio for large $|x|$; what is not as well-known is that the continued fraction is also useful for numerical evaluation.
Lentz, Thompson and Barnett derived an algorithm for numerically evaluating a continued fraction as an infinite product, which is more efficient than the usual approach of computing a continued fraction "backwards". Instead of displaying the general algorithm, I'll show how it specializes to the computation of Mills's ratio:
$\displaystyle Y_0=x,\,C_0=Y_0,\,D_0=0$
$\text{repeat for }j=1,2,\dots$
$$D_j=\frac1{x+jD_{j-1}}$$
$$C_j=x+\frac{j}{C_{j-1}}$$
$$H_j=C_j D_j$$
$$Y_j=H_j Y_{j-1}$$
$\text{until }|H_j-1| < \text{tol}$
$\displaystyle R(x)=\frac1{Y_j}$
where $\text{tol}$ determines the accuracy.
The CF is useful where the previously mentioned series starts to converge slowly; you will have to experiment with determining the appropriate "break point" to switch from the series to the CF in your computing environment. There is also the alternative of using an asymptotic series instead of the Laplacian CF, but my experience is that the Laplacian CF is good enough for most applications.

Finally, if you don't need to compute the (complementary) error function very accurately (i.e., to only a few significant digits), there are compact approximations due to Serge Winitzki. Here is one of them:
$$R(x)\approx \frac{\sqrt{2\pi}+x(\pi-2)}{2+x\sqrt{2\pi}+x^2(\pi-2)}$$
This approximation has a maximum relative error of $1.84\times 10^{-2}$ and becomes more accurate as $x$ increases.
