How can one calculate the probability of large random deviates (Z scores), such as 15 sigma from the expected mean?

Three references.  First, there is a nifty FORTRAN subroutine for the
complementary error function given in _Numerical_Recipes_ by Press,
Flannery, Teukolsky, and Vetterling (1986).  If we let G(t)=Pr(Z>t),
then the relationship between G() and erfc() is given by

    G(t)= (1/2) * erfc( t / sqrt(2) ),       t>0.

This approximation has a relative error of less than 10^-6, that is, it
is good to 6 significant figures.  Second, there is an excellent
algorithm due to Cooper (_Applied_Statistics_, vol 17, 186-187 [1968])
that is highly accurate.  Third, there is a great book called
_Applied_Statistics_Algorithms_ (1985), edited by Griffiths and Hill,
in which appears a short FORTRAN subroutine that produces 9-figure
accuracy (Hill, pp 126-129).

I have a book coming out late this year or early next entitled
_Elements_of_Statistical_Computing_, published by Chapman & Hall, in
which I talk about methods for constructing such approximations (the
references above are cribbed from my page proofs).  It also talks about
computing tail areas for other distributions than the normal, as well
as the inverse problem (finding a z value that corresponds to a given
probability).  One thing to note if you intend to "roll your own"
algorithm for this problem:  Out at z=15 you are much better off
computing the _logarithm_ of the probability rather than the
probability itself, since the probability is getting small at the rate
of about exp(-z^2/2).

-- Ron Thisted <thisted@galton.uchicago.edu>, Dept of Statistics, The
   University of Chicago

[ Editor: see also probability/t.distribution. ]