The computation of the multivariate normal distribution function is a common problem for statistical analysis in many different applications. The input to this problem is an m × m covariance matrix &Sgr; and the m lower and m upper integration limits for the m-dimensional integral of the multivariate normal density function. These integrals have been considered difficult to compute, particularly when m is large or high accuracy is required. Drezner’s paper describes a method for this problem when all lower integration limits are fixed. The method first converts the problem to a sum of simpler problems with all upper integration limits negative and then uses a sequence of appropriately chosen product Gauss integration rules to provide a sequence of successively more accurate approximations to these integrals. Results from a limited series of tests for the method are compared with results from Schervish’s MULNOR [1]. The new method is much faster for these tests, but still requires exponentially increasing computation time as m increases. While the new method appears to offer a significant improvement over Schervish’s method, two methods that are both probably much faster than either Drezner’s or Schervish’s methods are available. These other methods, developed by Deák [2] and me [3], will usually compute multivariate normal probabilities accurate to three or four decimal digits in times that appear to have only polynomial increase with m, for m as large as 15.