How to compute the probability associated with absurdly large Z-scores? Software packages for network motif detection can return enormously high Z-scores (the highest I've seen is 600,000+, but Z-scores of more than 100 are quite common).  I plan to show that these Z-scores are bogus.
Huge Z-scores correspond to extremely low associated probabilities.  The values of the associated probabilities are given on e.g. the normal distribution wikipedia page (and probably every stats textbook) for Z-scores of up to 6.  So...

Question: How does one compute the error function $1-\mathrm{erf}(n/\sqrt{2})$ for n up to 1,000,000, say?

I'm particularly after an already implemented package for this (if possible).  The best I've found so far is WolframAlpha, which manages to compute it for n=150 (here).
 A: You can approximate it with much simpler functions - see this Wikipedia section for more information.  The basic approximation is that $\textrm{erf}(x) \approx \textrm{sgn}(x)\sqrt{1 - \exp(-x^2 \frac{4/\pi + ax^2}{1+ax^2}})$
The article has an incorrect link for that section.  The PDF referenced can be found in Sergei Winitzki's files - or at this link.
A: A simple upper bound
For very large values of the argument in the calculation of upper tail
probability of a normal, excellent bounds exist that are probably as
good as one will get using any other methods with double-precision
floating point. For $z > 0$, let
$$
\renewcommand{\Pr}{\mathbb{P}}\newcommand{\rd}{\mathrm{d}}
S(z) := \Pr(Z > z) = \int_z^\infty \varphi(z) \rd z \>,
$$
where $\varphi(z) = (2\pi)^{-1/2} e^{-z^2/2}$ is the standard normal
pdf. I've used the notation $S(z)$ in deference to the standard
notation in survival analysis. In engineering contexts, they call this
function the $Q$-function and denote it by $Q(z)$.
Then, a very simple, elementary upper bound is
$$\newcommand{\Su}{\hat{S}_u} \newcommand{\Sl}{\hat{S}_\ell}
S(z) \leq \frac{\varphi(z)}{z} =: \Su(z) \> ,
$$
where the notation on the right-hand side indicates this is an
upper-bound estimate. This answer 
gives a proof of the bound.
There are several nice complementary lower bounds as well. One of the
handiest and easiest to derive is the bound
$$
S(z) \geq \frac{z}{z^2+1} \varphi(z) =: \Sl(z)  \> .
$$
There are at least three separate methods for deriving this bound. A
rough sketch of one such method can be found in this answer to a
related question.
A picture
Below is a plot of the two bounds (in grey) along with the actual
function $S(z)$.

How good is it?
From the plot, it seems that the bounds become quite tight even for
moderately large $z$. We might ask ourselves how tight they are and
what sort of quantitative statement in that regard can be made.
One useful measure of tightness is the absolute relative error
$$
\newcommand{\err}{\mathcal{E}} \err(z) = \left|\frac{\Su(z) -
S(z)}{S(z)}\right| \>.  
$$ 
This gives you the proportional error of
the estimate.
Now, note that, since all of the involved functions are nonnegative,
by using the bounding properties of $\Su(z)$ and $\Sl(z)$, we get
$$
\err(z) = \frac{\Su(z) - S(z)}{S(z)} \leq \frac{\Su(z) - \Sl(z)}{\Sl(z)} =
z^{-2} \> ,
$$
and so this provides a proof that
for $z \geq 10$ the upper-bound is correct to within 1%, for $z \geq 28$ it is correct to 
within 0.1% and for $z \geq 100$ it is correct to within 0.01%.
In fact, the simple form of the bounds provides a good check on other
"approximations". If, in the numerical calculation of more complicated
approximations, we get a value outside these bounds, we can simply
"correct" it to take the value of, e.g., the upper bound provided here.
There are many refinements of these bounds. The Laplace bounds mentioned here provide a nice sequence of upper and lower bounds on $S(z)$ of the form $R(z) \varphi(z)$ where $R(z)$ is a rational function.
Finally, here is another somewhat-related question and answer.
