The probability integral of the multivariate normal distribution (ND) has received considerable attention since W. F. Sheppard's (1900) and K. Pearson's (1901) seminal work on the bivariate ND. This paper evaluates the formula that represents the "n x n" correlation matrix of the "chi(sub i)" and the standardized multivariate normal density function. C. W. Dunnett and M. Sobel's formula for the univariate ND function, and R. E. Bohrer and M. J. Schervish's error-bounded algorithm for evaluating "F(sub n)" for general "rho(sub ij)" are discussed. Computationally, the latter algorithm is restricted to "n = 7"; even at "n = 7", it can take up to 24 hours for it to compute a single probability with 10(sup -3) accuracy on a computer than is capable of about 1-2 million scalar floating point operations/second. This report presents a fast and general approximation (APX) for rectangular regions of the multivariate ND function based on C. E. Clark's (1961) APX to the moments of the maximum of "n" jointly normal random variables. The performance of this APX compared to special cases in which the exact results are known and error-bounded reduction formulae show that the APX's accuracy is adequate for many practical applications where multivariate normal probabilities are required. The computational speed of the Clark APX is unparalleled. The error bound for the APX is about 10(sup -3) regardless of dimensionality, and accuracy increases with increases in "rho". The Clark algorithm provides a generalization of Dunnett's (1955) results to the case of general "rho(sub ij)", a natural application of which would be a generalization of Dunnett's test to the case of unequal sample sizes among the "k + 1" groups (i.e., multiple treatment groups compared to a single control group). One data table is included. (RLC)