{'fraction_non_alphanumeric': 0.030000656469506993, 'fraction_numerical': 0.011849274601194774, 'mean_word_length': 4.371473906911143, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
AMGKQ: An Efficient Implementation of Adaptive Multivariate Gauss-Kronrod Quadrature for Simultaneous Integrands in Octave/MATLAB

The algorithm AMGKQ for adaptive multivariate Gauss-Kronrod quadrature over hyper-rectangular regions of arbitrary dimensionality is proposed and implemented in Octave/MATLAB. It can approximate numerically any number of integrals over a common domain simultaneously. Improper integrals are addressed through singularity weakening coordinate transformations. Internal singularities are addressed through the use of breakpoints. Its accuracy performance is investigated thoroughly, and its running time is compared to other commonly available routines in two and three dimensions. Its running time can be several orders of magnitude faster than recursively called quadrature routines. Its performance is limited only by the memory structure of its operating environment. Included with the software are numerous examples of its invocation.


Introduction
As Press et al. [1992] state: "Integrals of functions of several variables, over regions with dimension greater than one, are not easy." We aim to make it so, at least for regions given by a hyper-rectangular volume in an arbitrary number of dimensions. Since one often requires the integration of several functions over the same region, an algorithm that performs the evaluations simultaneously is inherently more efficient than repeating the quadrature independently for each function. The use of a vectorized computing language, such as Octave or MATLAB, allows one to implement these evaluations with a minimum of coding.
In this article we propose AMGKQ, abbreviation for Adaptive Multivariate Gauss-Kronrod Quadrature. The main contributions of this work can be summarized as follows: • Vector, matrix, and binary singleton expansion operations are leveraged for efficiency.
• Variable transformations are used for improper integrals in multiple dimensions.
• An arbitrary number of breakpoints in multiple dimensions is possible.
The implementation AMGKQ.M is based on ADAPT.M by Alan Genz [van Dooren and de Ridder, 1976, Genz and Malik, 1980, Berntsen et al., 1991 and QUADGK.M by David Bateman [Eaton et al., 2009], and it makes use of subroutines provided by Walter Gautschi [1994Gautschi [ , 2004 and Randall LeVeque [2007]. The advice given by Shampine [2008] "to vectorize the evaluation of functions" is taken to extremes, as all the integrand values for each subregion are evaluated with a single call to the user supplied function. The quadrature coefficients are stored in persistent variables to reduce the computational load, and care is taken throughout to compute what is needed only once if possible. In the quest for efficiency, every FLOP counts.
The intended application of AMGKQ is in the context of Bayesian data analysis, where one usually finds a strongly peaked evidence density somewhere on the coordinate manifold against which expectation values of the observables are taken, but it should meet the requirements of a general purpose algorithm on par with those provided by Octave and MATLAB. Standard variable transformations are used when an improper integral is detected [Shampine, 2010], and internal singularities are avoided by use of breakpoints defining boundaries within the integration region. For good measure, it also does complex line (contour) integrals in the complex plane using the same machinery.

Statement of the Algorithm
In this section we will state the requirements of the algorithm, its initialization, and its main control loop. The basic theory of Gauss-Kronrod quadrature is assumed to be known to the reader, as are the implementations ADAPT.M and QUADGK.M. How the algorithm is used for contour integrals will be described at the end of this section.

Definition of the Integrand, the Region, and the Initial Subregions
The form of the user supplied function(s) F f (x, y, . . .) is important to the efficient implementation of its integration in a multivariate setting. Using the notation S for scalars, V for vectors, and M for matrices, what we require is Y = F(X), where X has size [N D , N X ] for N X points in N D dimensions and Y has size [N F , N X ] for N F integrands. Vectorized expressions and binary singleton expansion operators ⊕ and ⊗ should be used when coding F.
The region of integration is specified by the vectors A and B, each with N D elements, defining a hyperrectangular volume such that is the vector of N F results we are after with estimated error E. If necessary, the limits are swapped such that all A d < B d , accounting for any induced change of sign. Optionally, a matrix of breakpoints C with size [N D , N C ] can be supplied by the user; if none is given the default C is determined to be the midpoint of the region. The primary use of C is to inform AMGKQ of the locations of singular (or nearly singular) values of the integrand, but when doing complex line integrals C is used to define the path of the contour.
With C in hand, the initial subregions indexed by s are defined in terms of their central location L s and half-width H s as follows. That permutation of the N C locations in C which gives the shortest aggregate distance from A through C to B is selected. Starting near A, each point C c in turn is used to subdivide the region in which it is found into 2 N D partitions, discarding any null volumes. Optionally, the user can request that C be taken in the order originally specified. Points in C can be located on the outer boundary, C d,c = A d or C d,c = B d , or on internal boundaries without fubaring the initial subdivision. If all components in C are unique and not equal to any component of A or B, one has N C (2 N D − 1) + 1 subregions when finished, which sets the lower limit on the requested maximum number of subregions AMGKQ is allowed to consider. An example of the subdivision process in 2 dimensions is shown in Figure 1.

Gauss-Kronrod Quadrature in Multiple Dimensions
For any subregion labeled by s, we wish to compute the integral over the volume as efficiently as possible. To do so, we store in memory the abscissa in normalized units (between -1 and 1) for every contributing location X K of size [N D , N K ], as well as both their Gauss and Kronrod weights, W G and W K respectively. The (n G , n K ) Gauss-Kronrod quadrature rule pair (w G , w K ) in one dimension x K can be of any order n K = 2n G + 1, with tabulated values for those most commonly used [Holoborodko, 2011] and a double precision routine called for others [Laurie, 1997, Gautschi, 1994, 2004. In multiple dimensions N D > 1, one has N K = n N D K weights in W K and N G = n N D G weights in W G . Each element of W G and W K is the product of the weights in w G and w K respectively corresponding to the indexed location in X K , and X G is the subset of X K where every coordinate is of even parity. The abscissa locations in physical units X s can then be evaluated for the entire subregion by first calculating the locations along the central axes X L = (H s ⊗ x K ) ⊕ L s and then constructing X s from X L by indexing.
The integrand is then evaluated at all locations X s with a single call to the user supplied function, Y s = F(X s ). The Kronrod estimate of the integral is calculated with a matrix multiplication and a scalar multiplication d,s is the volume factor, and similarly for the Gauss The result for the subregion is given by the Kronrod estimate R s = Q K , and its variance is estimated as V s = (Q K − Q G ) 2 , taking the power along independent dimensions f . The values R s and V s are stored in memory for accumulation. The final result is the accumulation of all the subregion results R = ∑ s R s , and its estimated error is the square root of the accumulated subregion variances E = (∑ s V s ) 1/2 . Note that ADAPT multiplies its estimate of E by a factor of 3, while QUADGK accumulates the subregion errors V 1/2 s .

Selecting the Subregion and Direction for Subdivision
With each iteration of the main loop, that subregion among the N s present which has the single largest estimated variance across all integrands is selected for subdivision, while ADAPT selects the greatest error summed over integrands. When evaluating the direction for subdivision for N D > 1 according the magnitude of the fourth derivative of F, only that integrand with the largest variance is considered, while again ADAPT considers the sum over integrands. Neither method is particularly well-suited when the orders of magnitude of the integrands are vastly different, but how best to make relative the selection process for simultaneous integrands is not clear. When selecting the direction, AMGKQ focuses on the integrand that triggered the selection of the subregion.
The evaluation of the fourth derivative (in each dimension) is accomplished by using finite difference coefficients w 4 . Having evaluated Y s as part of the Gauss-Kronrod quadrature, those values are used again for this purpose. Since x K is not evenly spaced, the coefficients themselves must be calculated for the chosen order of quadrature rule [Fornberg, 1998, LeVeque, 2007. The abscissa locations X d along the central axes are identified, and a matrix multiplication yields the result F iv d,s = Y s (X d ) × w 4 in normalized units. That direction with the greatest fourth derivative in magnitude for the selected integrand is chosen for division by a factor of 2, such that each iteration contributes one additional subregion to the accumulation N s ← N s + 1. The calculation of R s , V s , and F iv d,s are implemented in the function gkint().

Convergence, Subregion Culling, and Termination Criteria
The user may request either or both an absolute tolerance E A and a relative tolerance E R for the convergence criterion. When all components of E are less than those of T, where T f is the greater of E A or E R |R f |, the algorithm considers itself globally converged. The estimated error E is a measure of the precision of the result, which is not quite the same thing as accuracy. To measure accuracy, one needs to know independently (analytically) what is the true value of the integral for comparison to its numeric approximation. One hopes, of course, that the precision and accuracy will be of the same order of magnitude, but difficult integrands can lead one to a result that is precisely wrong. The are two conditions under which a subregion may be culled from further consideration. The first is when its estimated error is sufficiently small as to not affect (the current estimate of) the final result. The second is when its half-width H s is approaching the limit of machine resolution in any dimension. Technically, that condition is met when subdividing a half-width would result in a subregion whose outermost Kronrod abscissa in physical units is indistinguishable numerically from the subregion boundary; otherwise, the integrand might be evaluated at the location of a breakpoint. When either condition is met, the contribution of the subregion to the accumulations is simply stored, and its location no longer considered. During testing, the second condition arose only when estimating ∞ 0 dx sin(x)/x = π/2, and the most accurate estimate of that integral was achieved by disabling the subregion culling entirely. If no subregions remain after culling, the algorithm considers itself converged but does send a unique flag to the user. The user also is warned if the second condition has been triggered.
There remain a few other conditions for which AMGKQ will terminate. When the maximum number of subregions N S requested by the user have been considered, the algorithm will return a flag along with its last values for the result and estimated error. Likewise, when a value of NaN or ±Inf is encountered, the algorithm will terminate with specific flags. These flags are meant to warn the user to inspect the result for accuracy. The user also is warned if any of the estimated errors do not meet the requested precision upon termination. The main loop of the algorithm can thus be stated as Algorithm 1.

ALGORITHM 1: Main Loop
Input: Initial subregions L s , H s , and function F(x). Output: Final estimates R, E, subregions evaluated N s , and flag I.
; also compute culling tolerances for H s ′′ and H s ′′′ ; end account for sign of R and take square root E = V 1/2 ; set flag I and express warnings;

Contour Integrals
Complex line (contour) integrals can be accomplished using the same machinery with no changes beyond some additional abs() functions that appear in the evaluation of the half-width tolerances. These integrals are restricted to N D = 1, which is understood to represent a single complex plane, and finite values for all components of A, B, and C. The points A and B are the starting and ending points of integration (which usually will be equal but are not required to be so), and the points in C determine the path of integration in a piecewise continuous linear fashion; no reordering of C is done in this case. Everything else proceeds the same as for the case of real integrals.

Improper Integrals and Variable Transformations
The algorithm AMGKQ can handle integrals that are improper, either because the integrand diverges at the boundary of the region or the domain of integration is itself unbounded. Singularities within the domain should be avoided by use of breakpoints. As long as the integrand is sufficiently well behaved, the result will be an accurate approximation to the value of the integral.

Edge Singularities
If a singular integrand is detected at either or both A and B, a variable transformation of the form dx F(x) = dy F(x y ) dx y /dy is employed to weaken the singularity; no attempt is made to weaken singularities at points in C. To find which dimensions are causing the singularity, a heuristic algorithm is employed. The idea is to define a point A ′ which is "near A", and similarly for B ′ . In case of infinite limits, one must check that "near A" is also "far from B" and act accordingly. Then, for each direction d, , and similarly for Y ′′ d,B . Any integrand which is not finite triggers the need for a variable transformation in that dimension at either or both endpoints as necessary. If the endpoints are themselves infinite in those dimensions, AMGKQ complains that the integral is divergent and throws an error. If not, the algorithm proceeds to effect the variable transformations. A schematic depiction for N D = 2 is shown in Figure 2.

Both A and B
This case is the most complicated thus will be considered first. Let {d ′ } ⊂ {d} be those dimensions for which singular integrands are detected at both A d ′ and B d ′ . The user has the option of selection either a trigonometric or a rational function for the transformation (in all d ′ , not independently). To accomplish the transformation, one needs to know x y ≡ x(y) and dx y /dy, as well as y x ≡ y(x) to find the new limits A d ′ and B d ′ and the breakpoints C d ′ in the new geometry. For the trigonometric transformation, those functions are . For the rational transformation we have where ; the imaginary part of y d ′ should be 0 to machine precision and can be discarded. For either transformation, one can construct the transformed integrand efficiently using binary singleton expansion such that F(

Just A or Just B
For these two cases only a rational transformation is available to the user. When {d ′′ } is the set of dimensions for which singular integrands are detected only at A d ′′ , the transformation functions are and when {d ′′′ } is the set of dimensions for which singular integrands are detected only at B d ′′′ , the transformation functions are All three cases are processed sequentially, which can lead to a final function of the form that gets passed to the main loop performing the actual quadrature, after accounting for any infinite limits that may be present.

Infinite Limits
When infinite limits appear in either A or B, a variable transformation is employed to map the manifold to a finite domain. The user has the option of selecting either a trigonometric or rational function for the transformation. Let {d ′′′′ } be the set of dimensions which have at least one infinite limit. For the trigonometric transformation, the required functions are and the required functions for the rational transformation are The infinite limit transformation F( is applied after any arising from edge singularities. If the default breakpoint C = (A + B)/2 is selected because none were specified by the user, one must account for infinite limits by taking C d ′′′′ = (A d ′′′′ + B d ′′′′ )/2 in the new coordinates.

Accuracy Testing
To test the accuracy of any numeric integration algorithm, one must assemble a collection of integrals whose values are known exactly. A well known collection is provided by John Burkardt 2009, 2011, a subset of which will be used here, sometimes modified for convenience. The analytic forms of these integrals have been included in the documentation provided with the code, as have example scripts that generate the results shown here. Some typos in the exact values quoted in the first library have been corrected, which are now evaluated in terms of their closed form solution.

Burkardt Tests
A set of 31 functions is selected for testing in one dimension as displayed in Table 1. The integrand is passed to AMGKQ in the form of an anonymous function handle. Default parameters of E A = eps(1) ≈ 1.5e-8, E R = 0, and N S = 2 N D × 100 control the algorithm, where eps(x) is the floating point resolution of value x. The limits for these functions are all finite. The number of subregions evaluated N s is displayed, as is the output flag: 2 means globally converged, 1 means locally converged, 0 means N S is reached, and <0 means Inf or NaN has been encountered. The estimated error is under the heading ERR, and the actual accuracy is under ACC.
For the vast majority of the selected functions, AMGKQ performs brilliantly. Only for the three functions numbered 29, 30, and 31 is ACC above the requested precision. The first is an example of a result that is precisely wrong, while the other two have an estimated error on par with their accuracy. The integration is repeated for these functions, which are displayed in Figure 3, with parameters E A = 0 and N S = 1000, and results are appended to the bottom of the table; what makes these integrands difficult are the sharp peaks in panel (a), the low-amplitude, high-frequency content in panel (b), and the wild oscillations in panel (c). Function number 29 is evaluated accurately when forced to converge locally, while function number 30 improves with more iterations. Function number 31 can be related to the sine integral which will be discussed later.
In two dimensions Burkardt [2011] does not give as many functions to investigate. These integrals all have limits which are the same for both directions, so only one value will be displayed for A and B. The results of the accuracy tests following the same procedure as above are shown in Table 2. Note that the two difficult integrands which were repeated both contain abs() as part of the operation; the discontinuity   (1/X) sin(1/X) 0.0 1.0 200 0 4.3e-01 7.2e-01 29 ∑ 3 k=1 (sech(10 k (X − k/5))) 2k 0.0 1.0 57 1 4.5e-17 0.0e+00 30 ∑ 40 k=1 cos(7 k Xπ/2)/2 k 0.0 1.0 1000 0 8.3e-05 -1.5e-05 31 (1/X) sin(1/X) 0.0 1.0 1000 0 1.4e-01 -2.2e-01 in the first derivative apparently makes high accuracy hard to achieve. Also note that the Gaussian function, which is representative of what is encountered in Bayesian data analysis, converges quickly compared to the others.
The numeric approximation of the integral of the sinc function f (x) = sin(x)/x over the semi-infinite domain is notoriously difficult. We have encountered one of its forms before among the difficult integrands of the Burkardt tests. If we ask QUADGK to approximate ∞ 0 dx sin(x)/x = π/2 ≈ 1.5708, it returns R = 5.7135 and E = 7.14403. If we call AMGKQ with N S = 1000 and disable subregion culling, we get a value of R = 1.5570, which has a relative accuracy of less than 1%. Speaking of relative accuracy, if we evaluate 15 10 dx sin(3x) cosh(x) sinh(x) ≈ 2.6e+10, with E A = 0 and E R = 1e-14, we find that the relative accuracy is indeed less than the requested relative precision.

Performance Testing
Having considered the advice given by Johnson [2002], let us embark on some performance testing. We will investigate the accuracy and running time in two and three dimensions of four functions, three of which are localized and one of which is oscillatory. In the order considered, the functions are a product over dimension of normal distributions , and a product of squared sinusoids F(x) = ∏ d sin 2 (x d ) cos 2 (x d ). The integration region begins as a square (or hyper-square) with sides of length 2 units, whose center is offset from the origin by up to half a unit in any direction, and is scaled by integer factors of k for successive runs. The normalization of each integral is set to unity for each k for consistency of comparison.
For two dimensional integrals N D = 2, Octave provides DBLQUAD, which calls recursively a chosen one dimensional quadrature routine such that vectorization of the integrand is only necessary for the first direction. For our purposes, we select QUADGK and QUADCC as our integration routines; QUADCC implements Clenshaw-Curtis quadrature rules. The adaptive Lobatto routine QUADL was abandoned for failing to terminate within a reasonable time when k becomes large, as was QUADV using an adaptive Simpson's rule. We also include ADAPT in our comparison to be complete, modified slightly to ignore the number of function evaluations and instead track N s . All quadrature routines are called with their default parameters for this test, except that they have E A = 1.0e-8 in common. The testing environment is Octave 3.8.1 using the ATLAS BLAS library running on a Pentium 4 CPU at 3.0 GHz with 3GB of RAM. The results are averaged over 3 trials to reduce their stochastic fluctuation.
We can see in Figure 4 that the accuracy of AMGKQ is comparable to that of DBLGK and DBLCC. Upon termination, all three produce a result whose accuracy is well below the requested precision. In contrast, ADAPT returns a result whose accuracy is on par with E A ; to be fair, that is all we asked for, but as we will see next, ADAPT has to work much harder to achieve a result that is nowhere near as accurate as the others. While Berntsen et al. [1991] recommend using a higher order rule for oscillatory integrands, the default rule of order 7 is used in ADAPT when evaluating the product of sinusoids in panel (d) for consistency of comparison. The relative performance of DBLGK and DBLCC varies with the choice of integrand, whereas AMGKQ is more consistent in that regard. Also shown is the accuracy of AMGKQ when it does all four integrands simultaneously.
In Figure 5 we compare the running times of the various implementations. Since the integrand functions are called in slightly different ways between DBLQUAD and AMGKQ, it is not really fair to compare the number of calls. The practical quantity which the user wants to minimize is running time, which is measured here in terms of CPU seconds, not wall time. For easy integrands (small k), AMGKQ returns a result up to   100 times faster than DBLGK, while for more difficult integrands the speedup factor is closer to 10, and it does not slow down as much as does DBLCC with increasing k. A factor of 10 might not seem like much for an operation that takes only a second, but in the context of Bayesian data analysis one often has to repeat variations of the same integral a large number of times. Furthermore, for real time predictive applications, every CPU cycle counts. We should mention that the initialization time for AMGKQ is not included in these comparisons. If it were, it would only affect k = 1 in panel (a) of Figure 5, raising its value by a factor of 10. After the first call, which need not produce anything useful, all the machinery for subsequent calls at the same order n G in the same number of dimensions N D is available in memory; since that feature is not part of the other implementations, it does not make sense to penalize AMGKQ for its inclusion. At any rate, all four integrands for each k can be evaluated by AMGKQ in one pass, thus doing them independently is already generous to its competition. The one pass running time displayed in the figure is one quarter of the time to do all four integrands simultaneously.
Let us next compare the number of subregions evaluated by AMGKQ and ADAPT, as shown in Figure 6. In panel (a) one sees that for the smallest k, AMGKQ converges after the initial subregions have been evaluated, while ADAPT requires on the order of 100 subregions (iterations) before it converges. For the larger k, there is a modest increase in N s upon termination for AMGKQ, while ADAPT reaches its maximum limit of N S = 2000 for the sinusoidal integrand. The vastly different values of N s for these two algorithms can only be explained by the superior performance of the Gauss-Kronrod quadrature rules.
We can repeat the comparison for N D = 3, at least for small values of k. Only a single set of runs is considered, on account of the length of time TRIPLEQUAD takes to converge. Likewise, ADAPT is no longer considered for the same reason. The accuracies shown in Figure 7 of TPLGK, TPLCC, and AMGKQ are all, as expected, well below the requested precision, as is the accuracy of the simultaneous integrands. What is interesting is the comparison of their running times, displayed in Figure 8. When the integrals are done independently, we see that AMGKQ outperforms TRIPLEQUAD by a factor greater than 100 and sometimes close to 1000. Interestingly, doing the integrals simultaneously appears to take slightly longer than their aggregate time, but the accuracy is not allowed to drift as much for large k. Obviously, performing multivariate quadrature with a recursive algorithm is not the quickest path to success.

Environment Limitations and the Curse of Dimensionality
As implemented, AMGKQ is not self-limiting; the number of dimensions N D , the number of integrands N F , and the order of quadrature rules (n G , n K ) can be arbitrarily large. Of course, there are practical limits im-    where N X is determined by the quadrature order (n G , n K ) and the number of dimensions N D . When variable transformations are in play, there can be a succession of functions that each create an array the size of X s in memory, and there must be room in RAM to hold them. Another limit is imposed by the class of indexing variable, which is implemented in Octave as a signed integer. On a 32-bit system, the maximum number of elements N 32 that can be stored in an array is one less than the maximum positive integer that can be represented, or N 32 = 2 31 − 2. If either product N D N X or N F N X is greater than N 32 , AMGKQ will fail to allocate room in memory for X s or Y s , respectively. On a fully 64-bit system (hardware and interpreter), the number of elements possible N 64 is much greater. The relation between the number of abscissae and the order of quadrature is N X = (2n G + 1) N D , such that a greater number of dimensions is feasible if one reduces the quadrature order.
There is simply no getting around the fact that the accurate numeric approximation of the integral of a multivariate function requires a lot of work. As a final test of AMGKQ, let us evaluate at order n G = 7 an integrand comprised of the product of the first few Burkardt tests in one dimension, with each additional function evaluated from an independent variable. As we can see in Table 4 and Figure 9, the number of elements in X s and Y s grows quite quickly with dimension. The running time shown in panel (c) includes the initialization time, since we are comparing apples to apples here, and demonstrates the curse of dimensionality. Not only must AMGKQ evaluate a strongly growing number of elements in X s and Y s as N D increases, but it also must do more work to evaluate each dimension's contribution to the integrand. Inspecting the column for N s in the table, we see that AMGKQ converged for all these integrals immediately after initialization. When the sixth function was appended, the Octave process consumed over 2GB of RAM and took more than a few minutes of wall time to evaluate. On a modern platform, the algorithm should be able to handle more dimensions than are considered here.

Outlook
During final preparations, the algorithm CHEBINT [Poppe and Cools, 2013] has come to our attention. This work focuses on the exposition of AMGKQ and its comparison to its parent algorithms ADAPT and QUADGK. It would be interesting in future work to compare the performance of AMGKQ and CHEBINT directly. Efficiency can be measured not only in terms of running time but also in terms of code complexity; AMGKQ accomplishes its goals with less than 1000 lines of code, including comments, examples, and nearly 200 lines of tabulated coefficients.
There remain opportunities to improve the efficiency of the implementation of AMGKQ. The most obvious upon reading the code is the manner in which the variable transformations are addressed. Rather than performing the transformations sequentially, it would be better to identify for each dimension the required transformation(s) and then effect the change of variable in a single function to reduce memory overhead and other expenses. It might also be better to transpose the sense in which X s and Y s are stored. Such detailed investigations of efficiency improvement are left for the interested reader to perform.

Conclusion
This work describes an efficient algorithm for the adaptive multivariate Gauss-Kronrod quadrature of simultaneous integrands and its implementation in Octave, AMGKQ. Its accuracy is comparable to the numerical integration routines provided by Octave, and its running time is much faster in multiple dimensions. Its performance is achieved by using vectorized code as much as possible, including in the user supplied integrand function. Its performance is limited only by the memory structure of its operating environment. The numeric approximation of integrals of functions of several variables might not be easy, but it has at least become easier.